{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SeqAcademy Multiomics Tutorial\n",
    "\n",
    "# Contents\n",
    "1. Installation\n",
    "    1. Set up channels\n",
    "    2. Create an environment and install the packages\n",
    "2. Alignment\n",
    "    1. HISAT \n",
    "    2. Samtools\n",
    "    3. MultiQC\n",
    "3. ChIP-Seq\n",
    "    1. MACS\n",
    "    2. Bedtools\n",
    "    3. IGV\n",
    "4. RNA-Seq\n",
    "    1. HTSeq\n",
    "    2. DESeq\n",
    "    3. Visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Installation\n",
    "\n",
    "Before running any programs, we'll make sure that each software is installed correctly. This tutorial uses Bioconda (https://bioconda.github.io/). Bioconda is a channel for the conda package manager specializing in bioinformatics software. The available packages are listed here: https://bioconda.github.io/recipes.html#recipes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set up channels\n",
    "\n",
    "You will need to add the bioconda channel as well as the other channels bioconda depends on. It is important to add them in this order so that the priority is set correctly (that is, bioconda is highest priority).\n",
    "\n",
    "The conda-forge channel contains many general-purpose packages not already found in the defaults channel. The r channel is only included due to backward compatibility. It is not mandatory, but without the r channel packages compiled against R 3.3.1 might not work."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tutorial uses cells written in python and unix to perform its analyses. Lines that are written in unix are prefixed by an exclamation point. \n",
    "\n",
    "Select the following cell and run it. To run a cell, select the cell, click \"Cell\" the upper taskbar, and select \"Run Cells\". Or click the cell and press shift + enter. Alternatively, the contents of any cell may be copy+pasted into the terminal emulator to run.\n",
    "\n",
    "The `--add channels` is an option that is supplied to the command ot tell it to add certain channels that you specify. The way it is written, the channels \"defaults\", \"conda-forge\", and \"bioconda\", would be added in that order. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!conda config --add channels defaults\n",
    "!conda config --add channels conda-forge\n",
    "!conda config --add channels bioconda"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create an environment and install the packages\n",
    "\n",
    "In this tutorial we will create an environment named \"tutorial\" and install the packages in there. Environments offer ways of installing packages in specific environments so they can be managed and run for different specifications. You can create, export, list, remove and update environments that have different versions of Python and/or packages installed in them. Switching or moving between environments is called activating the environment. You can also share an environment file.\n",
    "\n",
    "This command will create an environment \"tutorial\" in which to install the packages used in this tutorial.\n",
    "\n",
    "Run the following commands to create the environment. The `-n` flag specifies the name of the environment to create (which is called \"tutorial\") and the list of packages following the name are the packages that will be installed in the \"tutorial\" environment.\n",
    "\n",
    "This will most likely take 10-15 minutes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!conda create -n tutorial hisat2 multiqc macs2 bioconductor-deseq matplotlib ggplot samtools bioconductor-rsamtools bedtools htseq  --yes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then activate the environment with the following command."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# For Mac and Linux\n",
    "\n",
    "!source activate tutorial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# For Windows\n",
    "\n",
    "!activate tutorial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Alignment\n",
    "## HISAT (Hierarchical Indexing for Spliced Alignment of Transcripts)\n",
    "\n",
    "In this tutorial, we'll use Hisat to align the sample reads to a reference genome. Hisat automatically downloads and preprocesses the reads so they're ready to be aligned. Hisat (hierarchical indexing for spliced alignment of transcripts) is a highly efficient system for aligning reads from RNA sequencing experiments. HISAT uses an indexing scheme based on the Burrows-Wheeler transform and the Ferragina-Manzini (FM) index, employing two types of indexes for alignment: a whole-genome FM index to anchor each alignment and numerous local FM indexes for very rapid extensions of these alignments. HISAT’s hierarchical index for the human genome contains 48,000 local FM indexes, each representing a genomic region of ~64,000 bp."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The RNA-Seq data we'll use is from https://www.ncbi.nlm.nih.gov/Traces/study/?acc=SRP106028 and the ChIP-Seq data is from https://www.ncbi.nlm.nih.gov/Traces/study/?acc=SRP132584\n",
    "\n",
    "The model organism for this project is Yeast i.e. Saccharomyces cerevisiae. For RNA-Seq, yeast data between euploid and aneuoploid conditions will be compared. For ChIP-SEq, yeast data between 3AT-treated and untreated conditions will be compared.\n",
    "\n",
    "The following cell contains python code. \n",
    "\n",
    "The Python programming language comes with a variety of built-in functions. Among these are several common functions, including:\n",
    "\n",
    "+ print() which prints expressions out\n",
    "+ abs() which returns the absolute value of a number\n",
    "+ int() which converts another data type to an integer\n",
    "+ len() which returns the length of a sequence or collection\n",
    "\n",
    "These built-in functions, however, are limited, and we can make use of modules to make more sophisticated programs.\n",
    "\n",
    "Modules are Python .py files that consist of Python code. Any Python file can be referenced as a module. A Python file called hello.py has the module name of hello that can be imported into other Python files or used on the Python command line interpreter. \n",
    "\n",
    "Run the following cell to see what it does and observe the output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pandas import read_csv\n",
    "\n",
    "RNASeqSRARunTableFile='data/RNASeqSRA.tsv'\n",
    "RNASeqSRATable = read_csv(RNASeqSRARunTableFile, delimiter='\\t')\n",
    "RNASeqoutrun = (RNASeqSRATable[\"Run\"]).astype(list)\n",
    "RNASeqoutputSam = \"test/\" + RNASeqoutrun + \".sam\"\n",
    "RNASeqoutputAlignmentSummary = \"test/\" + RNASeqoutrun + \".txt\"\n",
    "RNASeqoutputMetrics = \"test/\" + RNASeqoutrun + \".metrics\"\n",
    "RNASeqoutputSortBam = \"test/\" + RNASeqoutrun + \".sorted.bam\"\n",
    "\n",
    "ChIPSeqSRARunTableFile='data/ChIPSeqSRA.tsv'\n",
    "ChIPSeqSRATable = read_csv(ChIPSeqSRARunTableFile, delimiter='\\t')\n",
    "ChIPSeqoutrun = (ChIPSeqSRATable[\"Run\"]).astype(list)\n",
    "ChIPSeqoutputSam = \"test/\" + ChIPSeqoutrun + \".sam\"\n",
    "ChIPSeqoutputAlignmentSummary = \"test/\" + ChIPSeqoutrun + \".txt\"\n",
    "ChIPSeqoutputMetrics = \"test/\" + ChIPSeqoutrun + \".metrics\"\n",
    "ChIPSeqoutputSortBam = \"test/\" + ChIPSeqoutrun + \".sorted.bam\"\n",
    "\n",
    "print(\"RNA-Seq run \" + RNASeqoutrun)\n",
    "print(\"ChIP-Seq run \" + ChIPSeqoutrun)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then run the following command to create the yeast index. The following command is a bash script `makeYeastIndex.sh` in the directory `scripts`. When run, the script downloads sequences for the latest Yeast release from Ensembl. By default, it builds and index for just the base files, since alignments to those sequences are the most useful.  To change which categories are built by this script, edit the CHRS_TO_INDEX variable in the `scripts/makeYeastIndex.sh` file. \n",
    "\n",
    "We will also use the `mkdir` command to make a directory `test` that will be used for storing hte output from this tutoirial. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!bash scripts/makeYeastIndex.sh\n",
    "!mkdir test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Align the RNA-Seq samples using Hisat.\n",
    "\n",
    "This step will most likely take several hours. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for index, individual in enumerate(RNASeqoutrun):\n",
    "    run = RNASeqoutrun[index]\n",
    "    summary = RNASeqoutputAlignmentSummary[index] \n",
    "    metrics = RNASeqoutputMetrics[index]\n",
    "    sam = RNASeqoutputSam[index]\n",
    "    bam = RNASeqoutputSortBam[index]\n",
    "    !hisat2 -x yeast_index/genome --sra-acc $run --new-summary --summary-file $summary --met-file $metrics -S $sam"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Samtools \n",
    "\n",
    "We'll use samtools to sort the output files and convert them to bam files.\n",
    "\n",
    "Sort the output files and convert them to bam files. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for index, individual in enumerate(RNASeqoutrun):\n",
    "    run = RNASeqoutrun[index]\n",
    "    summary = RNASeqoutputAlignmentSummary[index] \n",
    "    metrics = RNASeqoutputMetrics[index]\n",
    "    sam = RNASeqoutputSam[index]\n",
    "    bam = RNASeqoutputSortBam[index]\n",
    "    !samtools view -bSF4 $sam | samtools sort -o $bam\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Do the same thing for ChIP-Seq samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for index, individual in enumerate(ChIPSeqoutrun):\n",
    "    run = ChIPSeqoutrun[index]\n",
    "    summary = ChIPSeqoutputAlignmentSummary[index] \n",
    "    metrics = ChIPSeqoutputMetrics[index]\n",
    "    sam = ChIPSeqoutputSam[index]\n",
    "    bam = ChIPSeqoutputSortBam[index]\n",
    "    index = \"yeast_index/genome\"\n",
    "    !hisat2 -p 2 -x $index --sra-acc $run --new-summary --summary-file $summary --met-file $metrics -S $sam"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for index, individual in enumerate(ChIPSeqoutrun):\n",
    "    run = ChIPSeqoutrun[index]\n",
    "    summary = ChIPSeqoutputAlignmentSummary[index] \n",
    "    metrics = ChIPSeqoutputMetrics[index]\n",
    "    sam = ChIPSeqoutputSam[index]\n",
    "    bam = ChIPSeqoutputSortBam[index]\n",
    "    !samtools view -bSF4 $sam | samtools sort -o $bam"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MultiQC\n",
    "\n",
    "This section details quality control checks on the read data from either RNAseq or ChIPseq data using MultiQC. MultiQC takes all output and log files from an alignment software program and aggregates the information from all samples into one convenient report (html by default).\n",
    "\n",
    "MultiQC was installed earlier in the tutorial, so all we need to do is run it on the data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "MultiQC is configured to run the same no matter what type of sequencing data is available, therefore the same command can be used to analyze either our RNAseq data or our ChIPseq data.  We include the option 'hisat_output' since we are aligning using the HISAT2 program.  See http://multiqc.info/docs/ for more information."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use the 'hisat_output' option because we are analyzing data downloaded and aligned using the HISAT2 program.  We use the '--force' option to overwrite any previous versions of the multiqc_report.  '--quiet' only shows log warnings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!multiqc \"\".join(RNASeqoutrun) --quiet --outdir test/multiqc_rnaseq --force\n",
    "!multiqc \"\".join(ChIPSeqoutrun) --quiet --outdir test/multiqc_chipseq --force"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ChIP-Seq\n",
    "## MACS (Model-based Analysis for ChIP-Seq)\n",
    "\n",
    "Peak-calling is one of the main steps scientists use in determining the locations where protein is bound in DNA. Peak detection software, such as MACS (Model-Based Analysis for ChIP-Seq), call peaks using the aligned sequecnes as input and returns precise locations of predicted peaks as output. In this tutorial, we'll use MACS.\n",
    "\n",
    "More information about MACS: http://liulab.dfci.harvard.edu/MACS/Download.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ChIPSeqControl = ChIPSeqSRATable.loc[ChIPSeqSRATable[\"source_name\"] == \"Untreated\"][\"Run\"].astype(list)\n",
    "ChIPSeqTreatment = ChIPSeqSRATable.loc[ChIPSeqSRATable[\"source_name\"] != \"Untreated\"][\"Run\"].astype(list)\n",
    "\n",
    "print(\"Control sample \" + ChIPSeqControl)\n",
    "print(\"Treatment sample \" + ChIPSeqTreatment)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for index, individual in enumerate(ChIPSeqControl):\n",
    "    outputdirectory = \"test/\" + ChIPSeqTreatment.iloc[index]\n",
    "    name = ChIPSeqTreatment.iloc[index]\n",
    "    immunoprecipitate = \"test/\" + ChIPSeqTreatment.iloc[index] + \".sorted.bam\"\n",
    "    control = \"test/\" + ChIPSeqControl.iloc[index] + \".sorted.bam\"\n",
    "    !macs2 callpeak -c $control -t $immunoprecipitate -n $name --outdir $outputdirectory"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For an in-depth discussion of what MACS2 does: https://github.com/taoliu/MACS/wiki/Advanced:-Call-peaks-using-MACS2-subcommands"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bedtools\n",
    "\n",
    "In this tutorial, we'll use Bedtools to extract the intersecting regions of the MACS output between the experimental conditions.\n",
    "\n",
    "The Bedtools suite of programs is widely used for genomic interval manipulation or \"genome algebra\". \n",
    "\n",
    "First we'll sort the output. The following line uses the `sort` command to sort the MACS output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for index, individual in enumerate(ChIPSeqTreatment):\n",
    "    macs_output = \"test/\" + name + \"/\" + name + \"_peaks.narrowPeak\"\n",
    "    sort = \"test/\" + name + \"/\" + name + \"_peaks.narrowPeak.sorted\"\n",
    "    !sort -k 1,1 -k2,2n $macs_output > $sort"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we'll find the intersecting regions between the different experimental conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!bedtools intersect -a test/SRR6703661/SRR6703661_peaks.narrowPeak.sorted -b test/SRR6703663/SRR6703663_peaks.narrowPeak.sorted -u > test/ChIPSeqintersect.bed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Integrative Genomics Viewer\n",
    "\n",
    "A BAM file viewer will allow you to see your reads in an interactive graphical display. There are many different viewers available such as UCSC Genome Browser, Integrative Genomics Viewer (IGV), and NCBI Genome Workbench.\n",
    "\n",
    "To load .bam files please first do the following:\n",
    "\n",
    "Transfer the data to your computer. Amazon provides documentation on copying files with scp (Linux, OSX) and WinSCP (Windows). \n",
    "\n",
    "For scp, the command will be approximately: scp ubuntu@your-instance-name-or-ip-address:/home/ubuntu/data/transcript.fa\n",
    "\n",
    "Load annotation in IGV\n",
    "File -> Load from File ... -> ChIPSeqintersect.bed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# RNA-Seq\n",
    "## HTSeq (High-through sequencing)\n",
    "\n",
    "HTSeq is a Python library to facilitate the rapid development of RNA-Seq analysis. HTSeq offers parsers for many common data formats in HTS projects, as well as classes to represent data, such as genomic coordinates, sequences, sequencing reads, alignments, gene model information and variant calls, and provides data structures that allow for querying via genomic coordinates. In this tutorial we will use htseq-count, a tool developed with HTSeq that preprocesses RNA-Seq data for differential expression analysis by counting the overlap of reads with genes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pandas import read_csv\n",
    "\n",
    "gtf = \"data/Saccharomyces_cerevisiae.R64-1-1.84.gtf\"\n",
    "\n",
    "RNASeqSRARunTableFile='data/RNASeqSRA.tsv'\n",
    "RNASeqSRATable = read_csv(RNASeqSRARunTableFile, delimiter='\\t')\n",
    "RNASeqoutrun = (RNASeqSRATable[\"Run\"]).astype(list)\n",
    "RNASeqoutputSortBam = \"test/\" + RNASeqoutrun + \".sorted.bam\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for index, individual in enumerate(RNASeqoutputSortBam):\n",
    "    input = individual\n",
    "    output = individual + \".genecount.txt\"\n",
    "    !htseq-count -m intersection-nonempty -s no -f bam $input $gtf > $output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## DESeq (Differential Expression Sequencing)\n",
    "\n",
    "Estimate variance-mean dependence in count data from high-throughput sequencing assays and test for differential expression based on a model using the negative binomial distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!Rscript scripts/runDeseq.R"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualization\n",
    "\n",
    "The following script performs principal component analysis and creates volcano plots and bar graphs of RNA-Seq expression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!Rscript scripts/loadYeastGeneCounts.R"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.5.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}