module-2-chipseq/module-2-chipseq-d2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## CMM262: ChIP-seq workshop – Part 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Following the steps 1-5 on from Part 1, you should have a tagdirectory for each of the datasets (6 in total) and you should also have the UCSCfiles.   \n",
    "\n",
    "As we discussed, each dataset was created by a different antibody, and they can be divided into three types: TFs (transcription factors) HMs (histone modifications) and global input (background). 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.)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Run all homer commands in your terminal. First make sure to activate the r-bio conda environment. Make sure you are in the module-2-chipseq folder.**\n",
    "\n",
    "**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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "source activate r-bio\n",
    "cd ~/module-2-chipseq"
   ]
  },
  {
   "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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "mkdir chipseq_workshop/tagdirs/hms/\n",
    "mkdir chipseq_workshop/tagdirs/tfs/\n",
    "mkdir chipseq_workshop/tagdirs/input/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "mv chipseq_workshop/tagdirs/h3* chipseq_workshop/tagdirs/hms/\n",
    "mv chipseq_workshop/tagdirs/input-esc/ chipseq_workshop/tagdirs/input/\n",
    "mv chipseq_workshop/tagdirs/*esc chipseq_workshop/tagdirs/tfs/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "**6. One of the most common tasks with ChIP-seq data is to find ‘enriched’ regions commonly called “peaks”.** HOMER contains a command called findPeaks which is used to analyze tag directories for peaks. There are two common ways to use the command:\n",
    "\n",
    "```\n",
    "findPeaks <tag directory> -i <control tag directory> -style factor -o auto\n",
    "```\n",
    "\n",
    "or\n",
    "\n",
    "```\n",
    "findPeaks <path>/oct4-esc/ -i <path>/input-esc/ -style histone -o auto\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The difference between the two is in the “-style factor/histone” argument, which will tell the program to look for focal, fixed width peaks (factor; i.e. TF) vs. variable length peaks (histone; i.e. HM); the later is more common in the case of histone modifications. To find Oct4 peaks in the data, run the following command:\n",
    "\n",
    "```\n",
    "findPeaks <path>/oct4-esc/ -i <path>/input-esc/ -style factor -o auto\n",
    "```\n",
    "  \n",
    "This command will look for enriched regions and filter them based on several criterion, including ensuring that they have at least 4-fold more reads in peak regions relative to the control experiment (in this case `input-esc/`). The output will be stored in a HOMER-style peak file located in the Oct4 tag directory (`oct4-esc/peaks.txt`). The beginning of this file contains statistics and QC stats from the peak finding, including the number of peaks, number of peaks lost to input filtering, etc.     \n",
    "\n",
    "One field worth paying attention to is the **“Approximate IP efficiency”** which reports what fraction of reads from the experiment were actually found in peaks. For most decent experiments this value ranges from 1% to >30% (remember ChIP is an enrichment strategy... there is plenty of background in the data too!). Below this are the peaks along with enrichment statistics for each region.\n",
    "\n",
    "One other thing to note is that HOMER reports the results in a ‘peak’ file, which has a slightly different format from a traditional BED file format. To create a BED file from the peak file, use the tool pos2bed.pl (e.g. `pos2bed.pl oct4-esc/peaks.txt > oct4-esc.bed`. The “`> output.txt`” part at the end means that the results will be sent to stdout, and the “`> output.txt`” is used to capture the output information in a file.). BED files can be uploaded to IGV just like a bedGraph file. Also, most HOMER programs will work with either BED or peak files as input.\n",
    "\n",
    "\n",
    "Next we will find peaks for all samples using two ‘for loops’ – for the two types of data:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "for dir in chipseq_workshop/tagdirs/hms/*-esc; do\n",
    "    nohup findPeaks \"$dir\" -i chipseq_workshop/tagdirs/input/input-esc/ \\\n",
    "    -style histone -o auto 2> \"$dir\"/findPeaks.out &\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "for dir in chipseq_workshop/tagdirs/tfs/*-esc; do\n",
    "    nohup findPeaks \"$dir\" -i chipseq_workshop/tagdirs/input/input-esc/ \\\n",
    "    -style factor -o auto 2> \"$dir\"/findPeaks.out &\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make a directory for the annotation files (e.g., “annotations”) and convert the peak.txt and region.txt files to bed files:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "mkdir -p chipseq_workshop/annotations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "for dir in chipseq_workshop/tagdirs/hms/*-esc; do\n",
    "    dirname=\"$(basename $dir)\"\n",
    "    pos2bed.pl \\\n",
    "    \"$dir\"/regions.txt > chipseq_workshop/annotations/\"$dirname\".bed;\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "for dir in chipseq_workshop/tagdirs/tfs/*-esc; do\n",
    "    dirname=\"$(basename $dir)\";\n",
    "    pos2bed.pl \\\n",
    "    \"$dir\"/peaks.txt > chipseq_workshop/annotations/\"$dirname\".bed;\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:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "cp chipseq_workshop/annotations/* ~/private/viewer/annotations/"
   ]
  },
  {
   "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/annotations/*' ./\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",
    "Explore the original (previous to conversion to bed) peaks.txt and regions.txt files with less -S in the command line. For instance, each peak gets a score, see how the ones with high scores look on IGV vs. ones with low scores. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "**7. Now that we have identified peaks from our ChIP-seq data, it is time to figure out more information about where they are and what genes they might be regulating.** \n",
    "\n",
    "HOMER contains a program called annotatePeaks.pl that performs a wide variety of functions using peak/BED files. First, lets use it to perform basic annotation of the peak file. The annotatePeaks.pl program works like this:\n",
    "\n",
    "\n",
    "```\n",
    "annotatePeaks.pl <peak/BED file> <genome version> [options] > output.txt\n",
    "```\n",
    "\n",
    "To annotate peaks from the Oct4 experiment:\n",
    "\n",
    "```\n",
    "annotatePeaks.pl <path>/oct4-esc/peaks.txt mm9 > <path>/annotations/oct4.annotation.txt\n",
    "```\n",
    "\n",
    "You can do it using a for loop for the two types of datasets:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "for dir in chipseq_workshop/tagdirs/hms/*-esc; do\n",
    "    dirname=\"$(basename $dir)\"\n",
    "    nohup \\\n",
    "    sh -c \"annotatePeaks.pl ${dir}/regions.txt chipseq_workshop/mm9 > \\\n",
    "    chipseq_workshop/annotations/${dirname}.annotation.txt\" 2> \\\n",
    "    chipseq_workshop/annotations/annotatePeaks.${dirname}.out &\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and for tfs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "for dir in chipseq_workshop/tagdirs/tfs/*-esc; do\n",
    "    dirname=\"$(basename $dir)\"; \n",
    "    nohup \\\n",
    "    sh -c \"annotatePeaks.pl ${dir}/peaks.txt chipseq_workshop/mm9 > \\\n",
    "    chipseq_workshop/annotations/${dirname}.annotation.txt\" 2> \\\n",
    "    chipseq_workshop/annotations/annotatePeaks.${dirname}.out &\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we view the `“oct4.annotation.txt”` file with `less -S`, you’ll see several annotation columns. Take note of the columns specifying the nearest gene TSS, the distance, and the annotation of the genomic region the peak is located in. This annotation is split into two separate columns – one is basic (i.e., exon, promoter, intergenic, intron etc.), and a more detailed annotation that describes CpG islands, repeat elements, etc. You might have also noticed while the command was running that stats about annotation enrichment too."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "**8. The annotatePeaks.pl program can also be used to create histograms that display the read enrichment relative to given genomic features**, including transcription start sites (TSS) or any other set of regions the user wants to define. Since the TSS is so commonly used for this purpose, HOMER has a built-in annotation for TSS (based on RefSeq transcripts). The key parameters to create a histogram are the `-hist #` and `-size #` options, which control the binning size and total length of the histogram. The other important option is the `-d <tag directory>`, which specifies which experiments to compile histograms for. In general:\n",
    "\n",
    "```\n",
    "annotatePeaks.pl <peak/BED file> <genome version> -size <#> -hist <#> -d <Tag Directory> > output.txt\n",
    "```\n",
    "  \n",
    "(note that the peak/BED file can be replaced with the key word “tss” to make a histogram at the TSS). To create a histogram near the TSS with the data from the experiments we’ve looked at thus far, run the following:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "# should take ~30 secs\n",
    "annotatePeaks.pl tss chipseq_workshop/mm9 -size 8000 -hist 10 -d \\\n",
    "chipseq_workshop/tagdirs/hms/*-esc chipseq_workshop/tagdirs/tfs/*-esc \\\n",
    "chipseq_workshop/tagdirs/input/input-esc/ > \\\n",
    "chipseq_workshop/annotations/tss-histogram.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's check out the “tss-histogram.txt” file using pandas! You’ll notice that the first column gives the distance offsets from the TSS followed by columns corresponding to the ‘coverage’, ‘+ Tags’, and ‘- Tags’ for each experiment. Try graphing each as X-Y line graph using the first column as the X-coordinate to see the patterns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "histogram_df = pd.read_csv('~/module-2-chipseq/chipseq_workshop/annotations/tss-histogram.txt', sep='\\t', index_col=0)\n",
    "histogram_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "\n",
    "legend = []\n",
    "plt.figure(figsize=(15, 10))\n",
    "for column_name in histogram_df.columns:\n",
    "    if 'Coverage' in column_name:\n",
    "        histogram_df[column_name].plot()\n",
    "        plt.xlabel('Distance from TSS')\n",
    "        plt.ylabel('Read count')\n",
    "        legend.append(column_name)\n",
    "    plt.legend(legend)\n",
    "plt.title('Distance from TSS for each transcription factor or histone modification')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### Motif finding\n",
    "\n",
    "**9. DNA motif finding is a powerful technique to analyze ChIP-seq experiments.** Unlike gene expression data, ChIP-seq localizes signals to very specific regions of the genome allowing for accurate identification of the genetic signals responsible for recruiting various transcription factors. To use HOMER’s motif analysis program, run the findMotifsGenome.pl command using peak files from the experiments. In general the command works like this:\n",
    "\n",
    "```\n",
    "findMotifsGenome.pl <peak/BED file> <genome version> <output directory> [options]\n",
    "```\n",
    "\n",
    "Common options for motif finding are `-p <#cpu>` for parallel execution, and `-size <#>` to specify the size of the regions you which to search for DNA motifs. For transcription factor motifs, a good size is 100. Make a new directory (e.g., “motifs”). \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "mkdir -p chipseq_workshop/motifs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To find enriched motifs, copy the following for loop into your terminal to run it (again, make sure you are in the module-9-chipseq directory):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "# note that this will take ~7-8 mins to run\n",
    "for dir in chipseq_workshop/tagdirs/tfs/*; do\n",
    "    dirname=\"$(basename $dir)\"\n",
    "    mkdir -p chipseq_workshop/motifs/\"$dirname\"\n",
    "    nohup findMotifsGenome.pl \"$dir\"/peaks.txt \\\n",
    "    chipseq_workshop/mm9 chipseq_workshop/motifs/\"$dirname\" \\\n",
    "    -mask -size 100 -p 10 2> \\\n",
    "    chipseq_workshop/motifs/\"$dirname\"/findMotifsGenome.out &\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This command will perform several different steps, including checking for the enrichment of a library of known transcription factor motifs as well as perform a de novo search for enriched motifs. The “-mask” tells the program to mask repeat sequences in the genome. Depending on the speed of your server this program may take while to run. Once it’s finished, check out the homerResults.html files.\n",
    "\n",
    "This file will list the top de novo motifs found as well as provide stats and best matches to known transcription factor motifs. A second file called knownResults.html contains the enrichment statistics for known motifs. \n",
    "\n",
    "You might notice that both analyses indicate that Oct4 and Sox2 peaks in ESC are highly enriched for an OCT:SOX composite motif, which has been shown to be very important in establishing pluripotent enhancers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "**10. Now that you’ve found the most enriched motifs in a ChIP-seq experiment, it is worth it to see where those motifs are located (i.e., which peaks, etc.)**. One of the key outputs from motif finding are “motif files”, which contain the information needed to understand where the motif is located in the genome. For example, the top de novo motifs found during motif finding are located in the output directory in the homerResults/ directory. To make it easier, lets copy the top Oct4 motif to the file “topOct4.motif” in the main directory where we are executing these commands (alternatively you could save the motif file down from the HTML results):\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "mkdir -p chipseq_workshop/motifs/analysis\n",
    "cp chipseq_workshop/motifs/oct4-esc/homerResults/motif1.motif \\\n",
    "chipseq_workshop/motifs/analysis/topOct4.motif"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Now lets perform two separate analyses: first, lets create a histogram showing the motif positions relative to Oct4 peaks so we can check if it really looks enriched relative to the center of the Oct4 peaks. We can do this using the annotatePeaks.pl program like before, but instead of looking at read densities with the `-d <tag directory>` option, we can look at motif densities using the `-m <motif file>` option instead:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "annotatePeaks.pl chipseq_workshop/tagdirs/tfs/oct4-esc/peaks.txt \\\n",
    "chipseq_workshop/mm9/ -size 2000 -hist 10 \\\n",
    "-m chipseq_workshop/motifs/analysis/topOct4.motif > \\\n",
    "chipseq_workshop/motifs/analysis/output.Oct4motifHistogram.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's check out output.Oct4motifHistogram.txt, and create an X-Y plot to examine the distribution of the motif relative to the peaks.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "oct4motif_df = pd.read_csv('~/module-2-chipseq/chipseq_workshop/motifs/analysis/output.Oct4motifHistogram.txt', sep='\\t', index_col=0)\n",
    "oct4motif_df[oct4motif_df.columns[0]].plot()\n",
    "plt.xlabel(\"Distance from center\")\n",
    "plt.ylabel(\"Motif count\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, lets see where these motifs are located in the actual genome browser. To do this, we can run annotatePeaks.pl again, but this time we will not make a histogram and instead use the `-mbed <bedfile>` option to tell it to create a bed file of the motif positions that we can upload to the genome browser. Try the following:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "# will take ~30 secs\n",
    "annotatePeaks.pl chipseq_workshop/tagdirs/tfs/oct4-esc/peaks.txt \\\n",
    "chipseq_workshop/mm9/ -mbed chipseq_workshop/viewer/topOct4.motifTrack.bed \\\n",
    "-m chipseq_workshop/motifs/analysis/topOct4.motif > \\\n",
    "chipseq_workshop/viewer/peaks-topOct4motif.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now you can download and then load the `topOct4.motifTrack.bed` to IGV. In addition, the other output file `peaks-topOct4motif.txt` will contain the peak annotation results with an additional column showing the peaks that contain the given motif."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "cp chipseq_workshop/viewer/peaks-topOct4motif.txt ~/private/viewer/annotations/\n",
    "cp chipseq_workshop/viewer/topOct4.motifTrack.bed ~/private/viewer/annotations/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember to copy over the files onto your local computer using `scp`:\n",
    "```bash\n",
    "scp -r 'USERNAME@dsmlp-login.ucsd.edu:~/viewer/annotations/peaks-topOct4motif.txt' ./\n",
    "scp -r 'USERNAME@dsmlp-login.ucsd.edu:~/viewer/annotations/topOct4.motifTrack.bed' ./\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "**11. Finally, we want to gain some experience comparing ChIP-seq experiments.** At first pass it might make sense to make a Venn diagram comparing peaks from two experiments to see how many overlap. <ins>This is a horrible way to analyze ChIP-seq data!!!</ins> \n",
    "This is because many peaks are close to the threshold of detection, barely making the cut for statistical significance (or not) in one experiment or another. A good practice is to create a scatter plot comparing the read counts between two experiments directly at all of the sites where there is signal (i.e., peaks). To do this, first lets merge the peak files from the two experiments, collapsing peaks found that overlap:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "mkdir chipseq_workshop/peak-analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "mergePeaks chipseq_workshop/tagdirs/tfs/oct4-esc/peaks.txt \\\n",
    "chipseq_workshop/tagdirs/tfs/sox2-esc/peaks.txt > \\\n",
    "chipseq_workshop/peak-analysis/oct4andsox2.peaks.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When you run this command, you’ll notice that it will print out the numbers of overlapping and unique peaks from each file (i.e., Venn diagram). This can be useful to give you a general idea of how similar the experiments are, but be careful not to over interpret these values. Now that you have the combined features in the “oct4andsox2.peaks.txt” file, lets use annotatePeaks.pl to quantify the read counts from each experiment at each of the peaks by specify each tag directory with the -d option like the following:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "# will take ~30 secs\n",
    "annotatePeaks.pl chipseq_workshop/peak-analysis/oct4andsox2.peaks.txt \\\n",
    "chipseq_workshop/mm9 -d chipseq_workshop/tagdirs/tfs/oct4-esc/ \\\n",
    "chipseq_workshop/tagdirs/tfs/sox2-esc/ > \\\n",
    "chipseq_workshop/peak-analysis/oct4andsox2-scatter.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now open the oct4andsox2-scatter.txt file to directly look at the read counts (tag counts; last two columns) at each peak to find those with low levels in one experiment or the other. Try creating an X-Y scatter plot of the read counts with log-transformed axes to get a sense for how different (or similar) each experiment is from one another. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "scatter_df = pd.read_csv('~/module-2-chipseq/chipseq_workshop/peak-analysis/oct4andsox2-scatter.txt', sep='\\t')\n",
    "scatter_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "# Getting the last two columns...\n",
    "scatter_df.plot.scatter(x=scatter_df.columns[-2], y=scatter_df.columns[-1], edgecolors='white', size=40, linewidth=0.3)\n",
    "plt.xlabel('oct4-esc read counts')\n",
    "plt.ylabel('sox2-esc read counts')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have an appreciation for how similar the experiments are, lets try using the getDifferentialPeaks command to selective find peaks that are ‘specific’ to one experiment relative to another. In the following example we’ll specify “-F 2” to indicate that we want to find Sox2 peaks that are 2-fold higher in the Sox2 experiment relative to the Oct4 experiment:\n",
    "\n",
    "```\n",
    "getDifferentialPeaks <peak file to check> <target tag directory> <background tag directory> \n",
    "-F <fold change> > output.txt\n",
    "```\n",
    "\n",
    "and in our case:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "getDifferentialPeaks chipseq_workshop/tagdirs/tfs/sox2-esc/peaks.txt \\\n",
    "chipseq_workshop/tagdirs/tfs/sox2-esc \\\n",
    "chipseq_workshop/tagdirs/tfs/oct4-esc -F 2 > \\\n",
    "chipseq_workshop/peak-analysis/sox2-specific-peaks.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have found Sox2 peaks that are relatively uniquely bound by Sox2 (and not Oct4), look at the peak file using less -S, and look at several of the peaks in the genome browser to convince yourself that we have found interesting peaks.\n",
    "\n",
    "Finally, lets use motif finding to see if we can identify what is unique about the DNA sequences in the Sox2-specific peaks that might have lead to a differential recruitment of Sox2 versus Oct4 at these sites. Run motif finding on the Sox2 specific peaks:\n",
    "(reminder: this is the general command for finding peaks:\n",
    "\n",
    "```\n",
    "findMotifsGenome.pl <peak/BED file> <genome version> <output directory> [options]\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "mkdir -p chipseq_workshop/peak-analysis/motifs-sox2-specific"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "# will take ~3 mins\n",
    "nohup findMotifsGenome.pl \\\n",
    "chipseq_workshop/peak-analysis/sox2-specific-peaks.txt \\\n",
    "chipseq_workshop/mm9 chipseq_workshop/peak-analysis/motifs-sox2-specific \\\n",
    "-size 100 -p 10 -mask 2> \\\n",
    "chipseq_workshop/peak-analysis/findMotifsGenome-sox2.out &"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "**12. Notice anything different about these results relative to the results from all Sox2 peaks that we found in step 9?**"
   ]
  }
 ],
 "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}