module-2-chipseq/module-2-chipseq-d1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# BIOM262: ChIP-Seq workshop"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Introduction:**\n",
    "\n",
    "This workshop will walk you through an example of ChIP-seq analysis. We will focus on running tools as command lines and simple bash scripts – I recommend having a cheat-sheet like this [one](http://cheatsheetworld.com/programming/unix-linux-cheat-sheet/). \n",
    "\n",
    "\n",
    "We will use some common tools such as: \n",
    "* **bwa** for alignment (<a href=\"http://bio-bwa.sourceforge.net/bwa.shtml\" target=\"_blank\">http://bio-bwa.sourceforge.net/bwa.shtml</a>),\n",
    "* **IGV** for visualization (<a href=\"http://software.broadinstitute.org/software/igv/home\" target=\"_blank\">http://software.broadinstitute.org/software/igv/home</a>) \n",
    "* Most of the workshop will be done using **HOMER** (e.g., QC of the data, peak calling etc.; <a href=\"http://homer.ucsd.edu/homer/\" target=\"_blank\">http://homer.ucsd.edu/homer/</a>). HOMER was created by Chris Benner at UCSD, and I love the documentation and tutorials and the threaded humor. To install HOMER follow <a href=\"https://github.com/biom262/cmm262-2020/blob/master/Module_5/Notebooks/Install_Homer.ipynb\" target=\"_blank\">these instructions</a>."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note:** Please use the 16 CPU DataHub environment when running this notebook. Otherwise, some of the commands may take a while to run."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To make it smoother, first type the command below in your terminal window."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "source activate r-bio"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "During the workshop, and in general, it is always good to type the command and get the notes and use options of the command. Thus e.g., typing bwa would yield this output (capped after several lines):\n",
    "\n",
    "```\n",
    "(r-bio) agoren@dsmlp-jupyter-agoren$ bwa\n",
    "\n",
    "Program: bwa (alignment via Burrows-Wheeler transformation)\n",
    "Version: 0.7.17-r1188\n",
    "Contact: Heng Li <lh3@sanger.ac.uk>\n",
    "\n",
    "Usage:   bwa <command> [options]\n",
    "\n",
    "Command: index         index sequences in the FASTA format\n",
    "         mem           BWA-MEM algorithm\n",
    "         fastmap       identify super-maximal exact matches\n",
    "         pemerge       merge overlapping paired ends (EXPERIMENTAL)\n",
    "         aln           gapped/ungapped alignment\n",
    "         samse         generate alignment (single ended)\n",
    "         sampe         generate alignment (paired ended)\n",
    "         bwasw         BWA-SW for long queries\n",
    "\n",
    "         shm           manage indices in shared memory\n",
    "         fa2pac        convert FASTA to PAC format\n",
    "         pac2bwt       generate BWT from PAC\n",
    "         pac2bwtgen    alternative algorithm for generating BWT\n",
    "         bwtupdate     update .bwt to the new format\n",
    "         bwt2sa        generate SA from BWT and Occ\n",
    "\n",
    "Note: To use BWA, you need to first index the genome with `bwa index'.\n",
    "      There are three alignment algorithms in BWA: `mem', `bwasw', and\n",
    "      `aln/samse/sampe'. If you are not sure which to use, try `bwa mem'\n",
    "      first. Please `man ./bwa.1' for the manual.\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 0\n",
    "\n",
    "### 1. Organize directories\n",
    "\n",
    "Before we begin, we will create directories to organize our analysis:\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "# first, let's make sure we're in the right directory\n",
    "cd ~/module-2-chipseq\n",
    "mkdir chipseq_workshop\n",
    "mkdir chipseq_workshop/aligned\n",
    "mkdir chipseq_workshop/tagdirs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# persist the change in directory to other cells in this notebook\n",
    "%cd ~/module-2-chipseq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Generate symbolic links\n",
    "Let's generate a couple [symbolic links](https://linuxize.com/post/how-to-create-symbolic-links-in-linux-using-the-ln-command/) to make it easier to type file paths. You can think of them as shortcuts."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "ln -s path_to_file_or_directory symlink_name\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "ln -s ~/public/chipseq/fastqs/ chipseq_workshop/fastqs\n",
    "ln -s ~/public/chipseq/mm9/ chipseq_workshop/mm9\n",
    "ln -s ~/public/chipseq/mm9/mm9.fa.fai chipseq_workshop/fai"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Part I\n",
    "We will start with FASTQ files and perform many of the basic analysis tasks that one might normally do when analyzing ChIP-seq data. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **1.** Align FASTQ reads using `bwa mem`.\n",
    "The fastqs are at: `~/public/chipseq/fastqs`\n",
    "\n",
    "But we made a symbolic link so we can access them at `chipseq_workshop/fastqs`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "These files are originally from the following study investigating the roles that reprogramming factors play when transforming MEF (fibroblasts) into embryonic stem cells.\n",
    "[Chronis et al. Cooperative Binding of Transcription Factors Orchestrates Reprogramming](https://www.ncbi.nlm.nih.gov/pubmed/28111071)\n",
    "Sequencing Data: [GSE90893](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE90893)\n",
    "\n",
    "For this tutorial we extracted the ChIP-seq experiments for several transcription factors and histone modifications performed on ESC (embryonic stem cells). To reduce runtimes, only reads that mapped to chr17 (and chr17_random) are included. (the random chromosomes are explained in this link – http://genome.ucsc.edu/FAQ/FAQdownloads#download10.)\n",
    "\n",
    "**Question:** how would you generate such a file with only one chromosome?\n",
    "\n",
    "To align the reads we will use bwa. It is always a good practice to look at the manual of each tool you use, so you get an idea of options and documentation.\n",
    "Initially, to get an idea of the command we will run it on one file as follows:\n",
    "(notice if you have very long lines of code you can break it up with a \"\\\\\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "# this command will take ~33 secs to run\n",
    "bwa mem -t 16 chipseq_workshop/mm9/mm9.fa chipseq_workshop/fastqs/oct4-esc.chr17.2m.fastq \\\n",
    "| samtools sort -@ 16 -o chipseq_workshop/aligned/oct4-esc.chr17.2m.bam"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is a good practice to always double check datasets before you start analyzing them. For instance, use samtools to view the files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "samtools view chipseq_workshop/aligned/oct4-esc.chr17.2m.bam | head -n 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We should also validate that the files are indeed what they should be (e.g., aligned to chr 17, and have 2M reads). \n",
    "\n",
    "To calculate the number of reads do:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "samtools view chipseq_workshop/aligned/oct4-esc.chr17.2m.bam | wc -l"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, another option that is even better for here is to use samtools "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "samtools flagstat chipseq_workshop/aligned/oct4-esc.chr17.2m.bam"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you want to understand better the way SAM files are organized you can follow <a href=\"https://samtools.github.io/hts-specs/SAMv1.pdf\" target=\"_blank\">https://samtools.github.io/hts-specs/SAMv1.pdf</a> section 1.4.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "To do it properly, we will use a \"for loop\" to get BAMs from all the fastqs in the directory.\n",
    "\n",
    "Importantly, always check that what you think you asked to be done in your code, is indeed what the machine understood.\n",
    "echo is very useful - it should print out what the code is really saying. So first, let's make sure the files that are used in the for loop are the right ones:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "for f in chipseq_workshop/fastqs/*fastq; do\n",
    "    echo $f\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "note, if you want to view only a subset of the output (e.g., the directory has many files and you want to view only the first 3), you can just pipe the output to head or tail as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "for f in chipseq_workshop/fastqs/*fastq; do\n",
    "    echo $f\n",
    "done | head -n 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, let's get the file name for each file, and strip it from the path and anything else, such as the extension. We will also check it using echo:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "for f in chipseq_workshop/fastqs/*fastq; do\n",
    "    fname=\"$(basename $f .fastq)\"; \n",
    "    echo $fname\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's add the command we used before, but again, echo it as well. Given that the > symbol would write to the file, we will use quotes for the echo: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "for f in chipseq_workshop/fastqs/*fastq; do\n",
    "    fname=\"$(basename $f .fastq)\";\n",
    "    echo \"bwa mem -t 16 chipseq_workshop/mm9/mm9.fa $f | samtools sort -@ 16 -o > chipseq_workshop/aligned/$fname.bam\" ; \n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we can run it - need to make sure we removed the quotes and echo:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "# in total, this should take about ~2.7 mins to run\n",
    "for f in chipseq_workshop/fastqs/*fastq; do\n",
    "    fname=\"$(basename $f .fastq)\";\n",
    "    bwa mem -t 16 chipseq_workshop/mm9/mm9.fa $f | samtools sort -@ 16 -o chipseq_workshop/aligned/$fname.bam;\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This will produce BAM files for the 6 datasets. HOMER can analyze SAM files and if it receives BAM files it converts them to SAM so samtools has to be available (you can check that by typing “samtools” in the command line). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "--- \n",
    "### **2.** Create a “tag directory” \n",
    "\n",
    "**These commands should be run directly in the terminal**\n",
    "\n",
    "For the example Oct4 ChIP-seq experiment using the makeTagDirectory command. Start by typing makeTagDirectory (without any options) in your command line, it will provide the usage, some info about the command and a full list of program options – as I mentioned above, I highly recommend doing that whenever you use a new tool and a new command. \n",
    "\n",
    "Tag directories are analogous to sorted bam files and are the starting point for most HOMER operations like finding peaks, creating visualization files, or calculating read densities. The command also performs several quality control and parameter estimation calculations. The command has the following form:  \n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "makeTagDirectory <Output Tag Directory> [options] <input SAM/BAM file1> [input SAM/BAM file2] ...\n",
    "``` "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To create a tag directory for the Oct4 experiment, open the terminal, make sure you are in the module-9-chipseq folder, and run the following command with recommended options:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "# should take ~1 min\n",
    "makeTagDirectory chipseq_workshop/tagdirs/oct4-esc -genome chipseq_workshop/mm9 -checkGC \\\n",
    "chipseq_workshop/aligned/oct4-esc.chr17.2m.bam"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The command will take several seconds to run. What it is doing is parsing through the BAM file, removing reads that do not align to a unique position in the genome, separating reads by chromosome and sorting them by position, calculating how often reads appear in the same position to estimate the clonality (i.e. PCR duplication), calculating the relative distribution of reads relative to one another to estimate the ChIP-fragment length, calculating sequence properties and GC-content of the reads and performs a simple enrichment calculation to check if the experiment looks like a ChIP-seq experiment (vs. an RNA-seq experiment).\n",
    "\n",
    "The command creates a new directory, in this case named **oct4-esc**. Inside the directory are several text files that contain various QC results. \n",
    "\n",
    "Try looking at the following using the \"head\" command:\n",
    "\n",
    "> * **tagInfo.txt** - summary information from the experiment, including read totals.\n",
    "> * **tagFreqUniq.txt** - nucleotide frequencies relative to the 5’ end of the sequencing reads.\n",
    "> * **genomeGCcontent.txt** - distribution of ChIP-fragment GC%\n",
    "> * **tagAutocorrelation.txt** - relative distribution of reads found on the same strand vs. different strands.\n",
    "> * **tagCountDistribution.txt** - number of reads appearing at the same positions.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use this box to look at the various files\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "--- \n",
    "### **3.** Create “tag directories” for all samples\n",
    "\n",
    "By the following code, using a ‘for loop’, again pasting this command in the terminal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "# should take ~1 min\n",
    "for f in chipseq_workshop/aligned/*.bam; do\n",
    "    fname=\"$(basename $f .chr17.2m.bam)\";\n",
    "    nohup makeTagDirectory chipseq_workshop/tagdirs/\"$fname\" -genome chipseq_workshop/mm9 \\\n",
    "    -checkGC \"$f\" 2> chipseq_workshop/tagdirs/\"$fname\".out &\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this point you should have 6 tag directories. Look through the QC stats of the various ones."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "---\n",
    "### **4.** Next we will visualize the ChIP-seq experiments.\n",
    "\n",
    "By creating bedGraph files from the tag directories and using the IGV genome browser to look at the results. We will do this using the makeUCSCfile command. For most ChIP-seq experiments all you need to do is specify the tag directory and specify “-o auto” for the command to automatically save the bedGraph file inside the tag directory:\n",
    "\n",
    "```\n",
    "makeUCSCfile chipseq_workshop/tagdirs -o auto\n",
    "```\n",
    "\n",
    "For a specific dataset, e.g. Oct4, the command would be:\n",
    "\n",
    "```\n",
    "makeUCSCfile chipseq_workshop/tagdirs/oct4-esc/ -o auto\n",
    "```\n",
    "\n",
    "This creates the file “oct4-esc/oct4-esc.ucsc.bedGraph.gz”. This file format specifies the normalized read depth at variable intervals along the genome (use zmore and the filename to view the file format for yourself). \n",
    "\n",
    "Now make these for all samples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "mkdir -p chipseq_workshop/viewer\n",
    "for dir in chipseq_workshop/tagdirs/*-esc; do\n",
    "    dirname=\"$(basename $dir)\"\n",
    "    nohup makeUCSCfile $dir -o chipseq_workshop/viewer/$dirname.ucsc.bedGraph \\\n",
    "    2> chipseq_workshop/viewer/$dirname.out &\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### To view the file in the genome browser, do the following:\n",
    "\n",
    "Download the files to your computer. To do this, we will first need to copy the files to the private directory as follows:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "mkdir ~/private/viewer/\n",
    "cp chipseq_workshop/viewer/*gz ~/private/viewer/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, you can run the following command from the terminal of your local computer:\n",
    "```bash\n",
    "scp -r 'USERNAME@dsmlp-login.ucsd.edu:~/viewer/*gz' ./\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Note_: You must have `ssh` installed locally for this to work. Make sure to replace USERNAME with your UCSD AD username. Provide your AD password when prompted.\n",
    "\n",
    "**Open IGV.** Make sure you use the right genome (mm9; it is a good practice to have!) and drag the file to the center window (or select file -> load from file).\n",
    "\n",
    "The read pileups will display the relative density of ChIP-seq reads at each position in the genome. We only have data for chr17 in this example, so we can stick to that chromosome.\n",
    "\n",
    "---\n",
    "### **5.** See if there are any interesting patterns in the data that catch your eye.\n",
    "\n",
    "Try visiting the Pou5f1 locus (the gene for Oct4) by typing the gene name into the search bar at the top. Once at the Pou5f1 locus, zoom out (alt+click or scale on top right) to see if there any nearby sites that might resemble enhancers.\n",
    "\n",
    "Each dataset was created by a different antibody, and they can be divided into three types: TFs, HMs and global input. Since we will need to treat each type differently, I recommend making a directory for each – input, tfs and hms and move the tag directories to the relevant one (e.g. tfs/oct4-esc/, etc). \n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:r-bio]",
   "language": "python",
   "name": "conda-env-r-bio-py"
  },
  "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.8.12"
  },
  "nteract": {
   "version": "0.21.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}