How to assemble a transcriptome using NGS data.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 0 How to use this tutorial\n",
    "All hyperlinks are [clickable](img/IMG1.jpg \"What happened to the curious cat?\"), text in <span title=\"This tutorial was written in Jupyter (iPython notebook) and rendered on Github.\"><b>bold</b></span> has hovertext with additional information (move your mouse over the text to read the messages), and all code should be run in a Unix shell (Bash) terminal. Code will appear as indented blocks of text like this:\n",
    "```\n",
    "### These are comments within a block of code\n",
    "echo \"You can copy and paste this code in a Bash terminal on your computer\"\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1 Introduction\n",
    "<span title=\"The process of identifying the order of the individual components in a chain of molecules eg. nucleotides in RNA/DNA or amino acids in proteins.\"><b>Sequencing</b></span> plays an important role in diverse fields including \\[but not limited to\\] forensics, evolutionary biology, pharmacology, clinical pathology, oncology, and even anthropology. Nucleic acid sequencing is not exactly a *new* technology; scientists have been using various methods of nucleotide sequencing since the 1970s. Nearly half a centry later, the cost of sequencing has dropped dramatically thanks to numerous and discoveries and innovations.\n",
    "\n",
    "|<span title=\"Source: https://www.genome.gov/images/content/costpergenome2015_4.jpg\">The decreasing cost of sequencing</span>|<span title=\"Source: 'The cancer genome' by Michael R. Stratton, Peter J. Campbell, and P. Andrew Futreal. doi:10.1038/nature07943\">Evolution of sequencing technology</span>|\n",
    "|:-------:|:------:|\n",
    "|<span title=\"Source: https://www.genome.gov/images/content/costpergenome2015_4.jpg\"><img src='img/IMG2.jpg' width='400'></span>| <span title=\"Source: 'The cancer genome' by Michael R. Stratton, Peter J. Campbell, and P. Andrew Futreal. doi:10.1038/nature07943\"><img src=img/IMG4.jpg width='400'></span>|\n",
    "\n",
    "<center>[Click here](http://www.sciencedirect.com/science/article/pii/S0888754315300410) to learn more about the history sequencing technology.</center>\n",
    "\n",
    "Have you ever worked with sequencing data before? If your answer was '*Yes*', what were you trying to find? How was the data organized? Where did it come from? If your answer was '*No*', then let me welcome you to the era of !"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.1 Get to know your data\n",
    "\n",
    "In this tutorial we are going to be working with RNA sequencing data from the model organism *Saccharomyces cerevisiae*. As *S. cerevisiae* is a very well-studied organism, <span title=\"This isn't always a good idea; depending on your research objectives, you may need to consider which criteria were used to decide which features go in the reference annotations and what was discarded.\"><b>we will assume</b></span> that the reference genome & assembly are mostly correct and complete.\n",
    "\n",
    "What kind of organism is *S. cerevisiae*? *S. cerevisiae* is a budding yeast; a <span title=\"'Simple' is relative term...\"><b>simple</b></span> eukaryote with 16 chromosomes and <span title=\"According to yeastgenome.org, there are 5155 'validated' ORFs... but according to the NCBI, there are 6350 genes. Take note: annotations can vary between sources and releases.\"><b>~6000 genes</b></span>. Yeast are among the easiest organisms to grow and study in a laboratory setting. There are dozens of well-documented experimental techniques to manipulate the genomes of yeast- some of which you can even <span title=\"Try doing that with a chimpanzee... On second thought, please don't try that.\"><b>[do at home!](http://www.the-odin.com/diy-yeast-crispr-kit/)</b></span>\n",
    "\n",
    "|<span title=\"Cell cycle of budding yeast. Source: https://voer.edu.vn/file/3025\">Cell cycle</span>|<span title=\"S. cerevisiae dividing. Source: http://www.csb.ethz.ch/research/experimental-yeast-biology.html\">Your favorite sequence can be tagged with a fluorescent reporter!</span>|\n",
    "|:-------:|:-------:|\n",
    "|<span title=\"Cell cycle of budding yeast. Source: https://voer.edu.vn/file/3025\"><img src=\"img/IMG6.jpg\" width='350'>|<span title=\"S. cerevisiae dividing. Source: http://www.csb.ethz.ch/research/experimental-yeast-biology.html\"><img src=\"img/IMG5.png\" width='400'> </span>\n",
    "\n",
    "We grew our sample in a rich medium, meaning our yeast were happily growing before we extracted and purified their RNA. Our RNA sequencing (RNAseq) data is stored in a pair of FASTQ files: these files are composed of millions of reads which were run on the <span title=\"Additional information about our data: it is 50bp, paired-end, and stranded.\"><b>Illumina HiSeq2500 platform</b></span>-- a second generation sequencer.The main innovation of second generation sequencing [AKA next generation sequencing] was to split the task of identifying each base into a massively parallel process.<span title=\"Source: 'Next Generation Sequencing in Aquatic Models' by Yuan Lu, Yingjia Shen, Wesley Warren and Ronald Walter. DOI: 10.5772/61657\"><img src=img/IMG3.png width=500></span>\n",
    "If you haven't already, download the RNAseq data we will be working with, two FASTQ files called *s_cerevisiae_chrX_read1.fastq.gz* and * \ts_cerevisiae_chrX_read2.fastq.gz* from [this link](https://github.com/willblev/assembly_workshop_MA_2016). Typically a RNAseq run on a NGS platform will generate several gigabytes of data; our files are <span title=\"100MB vs 2.3GB\"><b> significantly smaller</b></span> because we are starting with only those reads which map to chromosome 10. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.2 Prerequesites\n",
    "\n",
    "First you need to download and <span title=\"You will need to download and install lots of different programs throughout your carrers as bioinformaticians. As each program author is different, the installation methods will vary. HINT: Most packages contain a file called README or INSTALL; this is always a good place to start.\"><b>install</b></span> the following programs:\n",
    "\n",
    "|   |   |\n",
    "|-------------|------------|\n",
    "|<li>[FastQC 0.11.5 for Linux](http://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.5.zip) or [FastQC 0.11.5 for OSX (Mac)](http://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.5.dmg)</li><br> <li> [Trimmomatic 0.36](http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.36.zip) (you need Java as well)</li> <br><li> [Bowtie2 2.2.9 for Linux](https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.2.9/bowtie2-2.2.9-linux-x86_64.zip/download) or [Bowtie2 2.2.9 for OSX (Mac)](https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.2.9/bowtie2-2.2.9-macos-x86_64.zip/download) </li><br><li> [SAMtools 1.3.1](https://github.com/samtools/samtools/releases/download/1.3.1/samtools-1.3.1.tar.bz2) </li><br><li>[Stringtie](https://ccb.jhu.edu/software/stringtie/#install) |<span title=\"XKCD: 'Will It Work' http://xkcd.com/1742/\"><img src=img/IMG9.png width=300></span>|\n",
    "\n",
    "After you have downloaded, decompressed, and installed all of these programs, you may need to tell your computer where to find the executeable files. When you try to execute a program in commandline, your terminal looks for executable files defined in the \\$PATH environment variable. To add the directory containing your new program to the \\$PATH environment variable, you can use the following code:\n",
    "```\n",
    "### for example, if the executeable file 'fastqc' is installed \n",
    "### in the directory '/home/username/Download/FastQC-0.11.5'\n",
    "\n",
    "export PATH=/home/username/Download/FastQC-0.11.5:$PATH\n",
    "```\n",
    "<br>\n",
    "You could also add this line to your ```~/.bash_profile``` if you wanted to store this path permanently. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# 2 Create a working directory and download reference files\n",
    "One of the most important things you will learn as a researcher is to organize and archive your work as you go, leaving notes and comments about your thought process and results. Use meaningful names for files using version control where necessary, and create directory structures that make sense for your data. As this tutorial is about RNA sequencing of *S. cerevisiae*, we will call our working directory \"S_cerevisiae_RNAseq_tutorial_2016\". Create this new directory on your computer wherever you please, then navigate to it.\n",
    "```\n",
    "### Make a new directory and navigate to it\n",
    "mkdir S_cerevisiae_RNAseq_tutorial_2016\n",
    "cd S_cerevisiae_RNAseq_tutorial_2016\n",
    "```\n",
    "Instead of using a browser to download the reference annotations like you did with with the programs you just installed, you can also download the files through a terminal window. \n",
    "```\n",
    "### Download the reference genome and annotations for Saccharomyces cerevisiae (S288c)\n",
    "wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/146/045/GCF_000146045.2_R64/GCF_000146045.2_R64_genomic.fna.gz \n",
    "wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/146/045/GCF_000146045.2_R64/GCF_000146045.2_R64_genomic.gff.gz \n",
    "```\n",
    "As these files are compressed, <span title=\"Some bioinformatic tools are designed to work with compressed files and thus do not require this step; as processors get faster and faster and data sets get larger and larger, working with compressed files will likely become the standard.\"><b>we have to decompress them</b></span>.\n",
    "```\n",
    "### Unzip the two files we just downloaded\n",
    "gunzip GCF_000146045.2_R64_genomic*.gz\n",
    "```\n",
    "Now let's rename these files so that they are easier to work with (and it will be easier to remember what they are). \n",
    "```\n",
    "### Rename the files with names are easier to remember\n",
    "mv GCF_000146045.2_R64_genomic.fna s_cerevisiae.fasta\n",
    "mv GCF_000146045.2_R64_genomic.gff s_cerevisiae.gff\n",
    "```\n",
    "We also need to download some other files to perform our analysis. These files are hosted on GitHub, so you can either clone the repository locally or download the files individually. \n",
    "```\n",
    "wget https://github.com/willblev/assembly_workshop_MA_2016/raw/master/s_cerevisiae_chrX_read1.fastq.gz\n",
    "wget https://github.com/willblev/assembly_workshop_MA_2016/raw/master/s_cerevisiae_chrX_read2.fastq.gz\n",
    "wget https://github.com/willblev/assembly_workshop_MA_2016/raw/master/s_cerevisiae_chrX.fasta\n",
    "wget https://github.com/willblev/assembly_workshop_MA_2016/raw/master/s_cerevisiae_chrX.gff\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.1 What did you just download? Take a closer look!\n",
    "Right now you should have a directory called 'S_cerevisiae_RNAseq_tutorial_2016' with 7 files in it. Make sure you have the same files with the following command:\n",
    "```\n",
    "### list the files in the directory \n",
    "ll -ht\n",
    "```\n",
    "You should see something like this:\n",
    "```\n",
    "total 114M\n",
    "-rw-r--r-- 1 username eg 459K Xxx xx xx:xx s_cerevisiae_chrX.gff\n",
    "-rw-r--r-- 1 username eg 738K Xxx xx xx:xx s_cerevisiae_chrX.fasta\n",
    "-rw-r--r-- 1 username eg  48M Xxx xx xx:xx s_cerevisiae_chrX_read2.fastq.gz\n",
    "-rw-r--r-- 1 username eg  46M Xxx xx xx:xx s_cerevisiae_chrX_read1.fastq.gz\n",
    "-rw-r--r-- 1 username eg 7.4M Xxx xx xx:xx s_cerevisiae.gff\n",
    "-rw-r--r-- 1 username eg  12M Xxx xx xx:xx s_cerevisiae.fasta\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 FASTA format\n",
    "The FASTA format is a common file format used for storing sequences; the more common extensions are: .fasta .fa .fna  (less common are .faa .frn and .ffn). \n",
    "It consists of a single line of identifying information, followed by several lines containing the sequence itself; one FASTA file can contain multiple sequences. \n",
    "The Fasta file *s_cerevisiae.fasta* contains the sequences of all 16 chromosomes. To see what this looks like, type:\n",
    "```\n",
    "head s_cerevisiae.fasta\n",
    "```\n",
    "You can see that the 'header' ie. the line containing the information about the sequence begins with the character *>* . We can see how many FASTA sequences are contained in this file by counting the number of lines which begin with *>*:\n",
    "```\n",
    "grep '^>' s_cerevisiae.fasta | wc -l\n",
    "```\n",
    "*Grep* looks through the file line-by-line, and returns all lines which match our regular expression '^>' (lines beginning with the character *>*). *Wc -l* counts the number of lines in the output from *grep*. Is this number different from the number of chromosomes *S. cerevisiae* is supposed to have? <span title='Usually when they ask this, you know that the answer is \"yes\"'><b>If so</b></span>, can you figure out why? Here's a <span title='If you use grep without piping the results to wc, you will see the content of each header. Is one of them a little different than the others?'><b>hint</b></span> if you need it.\n",
    "\n",
    "You may also notice that after *>* each header begins with NC_00xxxx.x- this is a [RefSeq identifier](http://www.ncbi.nlm.nih.gov/refseq/). It is a unique identifier, which helps to prevent mixing up 'ChrI' of *H. sapiens* and 'ChrI' of *D. melanogaster* if you are working with multiple organisms. \n",
    "\n",
    "Here is a table which translates the RefSeq IDs in our files to the common chromosome names:\n",
    "\n",
    "| ID in s_cerevisiae.fasta | Chromosome |\n",
    "|---------------------|------------|\n",
    "| NC_001133.9         | I          |\n",
    "| NC_001134.8         | II         |\n",
    "| NC_001135.5         | III        |\n",
    "| NC_001136.10        | IV         |\n",
    "| NC_001137.3         | V          |\n",
    "| NC_001138.5         | VI         |\n",
    "| NC_001139.9         | VII        |\n",
    "| NC_001140.6         | VIII       |\n",
    "| NC_001141.2         | IX         |\n",
    "| NC_001142.9         | X          |\n",
    "| NC_001143.9         | XI         |\n",
    "| NC_001144.5         | XII        |\n",
    "| NC_001145.3         | XIII       |\n",
    "| NC_001146.8         | XIV        |\n",
    "| NC_001147.6         | XV         |\n",
    "| NC_001148.4         | XVI        |\n",
    "| NC_001224.1         | Mito       |\n",
    "\n",
    "We have another FASTA file called *s_cerevisiae_chrX.fasta*; unsurprisingly, this file contains only the sequence of chromosome 10. You can check to make sure this FASTA file only has one sequence using *grep* and *wc* again:\n",
    "```\n",
    "grep '^>' s_cerevisiae_chrX.fasta | wc -l\n",
    "```\n",
    "\n",
    "Does the ReqSeq identifier in this file match the identifier for chromosome 10?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "### 2.3 GFF (General Feature Format)\n",
    "Another ubiquitous file format is the GFF (General Feature Format), which is often used to store information about annotated genomic features like genes, exons, ORFs, etc. Each line in a GFF is separated into 9 columns by tabs:\n",
    "\n",
    "1. seqname - name of the chromosome or scaffold; chromosome names can be given with or without the 'chr' prefix. Important note: the seqname must be one used within Ensembl, i.e. a standard chromosome name or an Ensembl identifier such as a scaffold ID, without any additional content such as species or assembly. See the example GFF output below.\n",
    "2. source - name of the program that generated this feature, or the data source (database or project name)\n",
    "3. feature - feature type name, e.g. Gene, Variation, Similarity\n",
    "4. start - Start position of the feature, with sequence numbering starting at 1.\n",
    "5. end - End position of the feature, with sequence numbering starting at 1.\n",
    "6. score - A floating point value.\n",
    "7. strand - defined as + (forward) or - (reverse).\n",
    "8. frame - One of '0', '1' or '2'. '0' indicates that the first base of the feature is the first base of a codon, '1' that the second base is the first base of a codon, and so on..\n",
    "9. attribute - A semicolon-separated list of tag-value pairs, providing additional information about each feature.\n",
    "\n",
    "Let's have a look at the reference annotations we downloaded from the NCBI:\n",
    "```\n",
    "head s_cerevisiae.gff\n",
    "```\n",
    "The first few lines of this file are commented out with *#* and are there to provide information about the source of the annotations. After these, each line describes a different feature, starting with the chromosome (see table above).  To see which features are included in our GFF, we can look more closely at the third field (feature type):\n",
    "```\n",
    "grep -v '^#' s_cerevisiae.gff |cut -f3 | sort | uniq -c\n",
    "```\n",
    "First we use the command *grep -v '^#'* to ignore all lines beginning with '#' ie. lines that are commented out, then we pipe the results to the program *cut* which selects only the third field, then we use *sort* to arrange the array alphabetically, and finally we use the program *uniq -c* to remove duplicate entries and count the number of times we saw each type of feature. \n",
    "\n",
    "How many coding sequences are there? How many genes? How many exons? Can you figure out how many genes are multi-exonic?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.4 FASTQ format\n",
    "The FASTQ format is typically used to record the results of a sequencing run. These runs produce millions of reads with a prediction of each base in the read and the quality of each prediction. Each read in a FASTQ file is written as four lines including a unique identifier, a sequence, and a per-base quality score:\n",
    "```\n",
    "1. Unique identifier        @HWI-D00021:29:C6KF5TBX:9:1821:1914:3311 2:N:0:GGCTAC\n",
    "2. Sequence                 TAGCTCTGGCACCAATTGGTTAACTCTTCTCCAAAGCAATGGGGTCAAAA\n",
    "3.                          +\n",
    "4. Quality score            BBBGGGGGGGGGGGGEBBGGGGGGBGGGGGFEGGGGGGGGGGGEGGGGGG\n",
    "\n",
    "```\n",
    "Have a look at the first 9 lines in one of our FASTQ files with *zcat* (like *cat*, but can be used to view zipped files) and *head*:\n",
    "```\n",
    "zcat s_cerevisiae_chrX_read1.fastq.gz | head -n 9\n",
    "```\n",
    "You may have noticed that we have *two* FASTQ files. This is because our RNAseq data is paired-end; paired-end is more expensive than single-end, but paired-end reads generally align better to your reference genome. Paired-end data is also more suitable than single-end data if you want to detect novel transcripts, new splicing isoforms, genomic rearrangements (such as inversions, deletions, or gene fusion events), and repetitive sequence elements. \n",
    "\n",
    "<span title=\"Source: http://www.yourgenome.org/sites/default/files/images/illustrations/bioinformatics_single-end_pair-end_reads_yourgenome.png\"><img src=\"img/IMG8.png\" width='560'></span> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3 Quality control with *FastQC*\n",
    "FastQC aims to provide a simple way to do some quality control checks on raw sequence data coming from high throughput sequencing pipelines. It provides a modular set of analyses which you can use to give a quick impression of whether your data has any problems of which you should be aware before doing any further analysis. It generates a HTML file with graphs that allow you to quickly spot artifacts or other problems with your raw data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Sometimes after downloading a program you may need to  <span title=\"Whenever you get a warning that says 'Permission denied' you know that the problem is with file permissions.\"><b>modify the permissions</b></span> of a file to be able to execute the program on your system. [Learn more about Linux file permissions here.](https://www.linux.com/learn/understanding-linux-file-permissions) We need to give the file \"fastqc\" permission to run on your system; you can give this file *executable* permission with the following command:\n",
    "```\n",
    "chmod a+x fastqc\n",
    "```\n",
    "Again, we first create a directory to contain the output of fastqc\n",
    "```\n",
    "### We create a new directory which will contain the results of our analysis\n",
    "mkdir fastqc\n",
    "```\n",
    "And then run the analysis\n",
    "```\n",
    "### Run fastqc on the raw sequencing reads [using 3 threads to speed things up]\n",
    "fastqc  s_cerevisiae_chrX_read1.fastq.gz -o fastqc -t 3 \n",
    "fastqc  s_cerevisiae_chrX_read2.fastq.gz -o fastqc -t 3 \n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 Analysis of *FastQC* results\n",
    "Open the files s_cerevisiae_chrX_read1_paired_fastqc.html and s_cerevisiae_chrX_read2_paired_fastqc.html in your browser. Does our data have any warnings or failures? Is the data <span title=\"In my experience, it is rare to have a perfect score with your raw data; that's why these checks are important! Depending on your project/objectives, warnings and errors from FastQC do not always require that your data be 'fixed'.\"><b>ready for analysis</b></span>? You can see [the description of each FastQC test and possible causes for failures here](http://www.bioinformatics.babraham.ac.uk/projects/fastqc/Help/3%20Analysis%20Modules/). The beginning of the FastQC results should look like this: <span title=\"FastQC results for s_cerevisiae_chrX_read1.fastq.gz  Per-base sequence quality\"><img src=\"img/IMG10.png\" width='560'></span> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4 Trimming RNAseq reads with *Trimmomatic*\n",
    "NGS sequencing typically results in <span title=\"Our samples had approximately 25 million read pairs.\"><b>tens of millions of reads</b></span>. A proportion of these reads will contain artifacts or low-quality bases which we would like to remove before starting our analyses. *Trimmomatic* performs a variety of useful trimming tasks for illumina paired-end and single ended data. The selection of trimming steps and their associated parameters are supplied on the command line. Our sample of *S. cerevisiae* was prepared with the TruSeq RNA library preparation protocol using adapter number 11, we will search through our reads (the FASTQ files) and remove this adapter. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Instead of downloading a file and editing it, we can simply create a file containing the one adapter sequence by running the follwing code:\n",
    "```\n",
    "### This code will create a file called illumina_adapter.fasta which contains\n",
    "### the sequence of TruSeq adapter number 11\n",
    "cat <<EOF > illumina_adapter.fasta\n",
    ">TruSeq_Adapter_Index_11\n",
    "GATCGGAAGAGCACACGTCTGAACTCCAGTCACGGCTACATCTCGTATGCCGTCTTCTGCTTG\n",
    "EOF\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "We need to create a directory to store the output of *Trimmomatic*\n",
    "```\n",
    "### First we will create a new directory within our working directory to store the trimmed fastq files\n",
    "mkdir trimmomatic\n",
    "```\n",
    "We can then run *Trimmomatic* on our RNAseq raw data. We use several parameters to tell the program how we want to filter and trim our data:\n",
    "* ILLUMINACLIP: Cut adapter and other illumina-specific sequences from the read.\n",
    "* LEADING: Cut bases off the start of a read, if below a threshold quality\n",
    "* TRAILING: Cut bases off the end of a read, if below a threshold quality\n",
    "* SLIDINGWINDOW: Perform a sliding window trimming, cutting once the average quality within the window falls below a threshold.\n",
    "* MINLEN: Drop the read if it is below a specified length\n",
    "\n",
    "```\n",
    "### Running Trimmomatic on our RNAseq data. You may need to add the directory/path to \n",
    "### the command if Java can't locate it\n",
    "java -jar Trimmomatic-0.36/trimmomatic-0.36.jar \\\n",
    "PE -phred33 s_cerevisiae_chrX_read1.fastq.gz s_cerevisiae_chrX_read2.fastq.gz \\\n",
    "trimmomatic/s_cerevisiae_chrX_read1_paired.fastq trimmomatic/s_cerevisiae_chrX_read1_unpaired.fastq \\\n",
    "trimmomatic/s_cerevisiae_chrX_read2_paired.fastq trimmomatic/s_cerevisiae_chrX_read2_unpaired.fastq \\\n",
    "ILLUMINACLIP:illumina_adapter.fasta:2:33:20:2:true LEADING:36 TRAILING:32 SLIDINGWINDOW:4:30 MINLEN:35\n",
    "```\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "\n",
    "We can then run the FastQC analysis again to see if we have fixed any of the warnings/failures we saw previously.\n",
    "```\n",
    "### Run fastqc on the pair of fastq files that trimmomatic produced [using 3 threads to speed things up]\n",
    "fastqc trimmomatic/s_cerevisiae_chrX_read1_paired.fastq -o fastqc -t 3 \n",
    "fastqc trimmomatic/s_cerevisiae_chrX_read2_paired.fastq -o fastqc -t 3 \n",
    "```\n",
    "\n",
    "**If you had problems running this section using trimmomatic, you can download the trimmed fastq files here:**\n",
    "* [s_cerevisiae_chrX_read1_paired.fastq.gz](https://github.com/willblev/assembly_workshop_MA_2016/raw/master/trimmomatic/s_cerevisiae_chrX_read1_paired.fastq.gz)\n",
    "* [s_cerevisiae_chrX_read2_paired.fastq.gz](https://github.com/willblev/assembly_workshop_MA_2016/raw/master/trimmomatic/s_cerevisiae_chrX_read2_paired.fastq.gz)\n",
    "\n",
    "You may need to move these files inside of your project folder and then into a directory called 'trimmomatic', as the next steps assume the files are in the proper location. You may also need to unzip these files with the command *gunzip*. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5 Map RNAseq reads to the reference genome with *Bowtie2*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "After we check to make sure that our raw RNAseq reads are useable, we can map them to our reference genome. In our case, we will be mapping only to a smaller part of the genome, chromosome X. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.1 Build an index\n",
    "First we must create a* bowtie2* index using chromosome 10 so that the millions of reads can be mapped efficiently to the reference sequence. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "```\n",
    "### Building a bowtie2 index with our .fasta file. It will create several files with the *.bt2 extension\n",
    "bowtie2-build  s_cerevisiae_chrX.fasta  s_cerevisiae_chrX\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.2 Map the reads with *Bowtie2*\n",
    "Now that we have built the index, we need to tell* bowtie2* where to find these index files:\n",
    "```\n",
    "### create env variable to show bowtie2 know where the index files are\n",
    "export BOWTIE2_INDEXES='.'\n",
    "```\n",
    "Then we can run *bowtie2 *to map our reads to the reference, which will result in a .SAM file\n",
    "```\n",
    "### Run bowtie2 to map our reads to the reference genome\n",
    "bowtie2 -x s_cerevisiae_chrX -p 2 \\\n",
    "-1 trimmomatic/s_cerevisiae_chrX_read1_paired.fastq \\\n",
    "-2 trimmomatic/s_cerevisiae_chrX_read2_paired.fastq \\\n",
    "-S 's_cerevisiae_chrX.sam'\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "To be able to use our mapped reads with *Stringtie*, we need to sort the mapped reads by location- we use the *samtools* suite of tools to convert the .SAM file into a .BAM file and sort it. A BAM file contains the all the same information as a SAM file, but it is <span title=\"BAM files are compressed using the BGZF format; they can be indexed to drastically improve the efficiency of certain operations and reduce overall file size.\"><b>compressed.</b></span> \n",
    "```\n",
    "### convert the .sam file into a .bam file and sort it \n",
    "samtools view -Su s_cerevisiae_chrX.sam > s_cerevisiae_chrX.bam\n",
    "samtools sort s_cerevisiae_chrX.bam  s_cerevisiae_chrX_sorted.bam\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 6. Assemble the mapped reads into a transcriptome with *Stringtie*\n",
    "*Stringtie* uses aligned RNA-Seq reads and assembles the alignments into a parsimonious set of transcripts. It uses a novel network flow algorithm as well as an optional de novo assembly step to assemble and quantitate full-length transcripts representing multiple splice variants for each gene locus.  <span title=\"Source: 'StringTie enables improved reconstruction of a transcriptome from RNA-seq reads' by Mihaela Pertea, Geo M Pertea, Corina M Antonescu, Tsung-Cheng Chang, Joshua T Mendell, and Steven L Salzberg. doi:10.1038/nbt.3122\"><img src='img/IMG7.jpg' width='750'></span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Create a directory to store the files generated by *Stringtie*\n",
    "```\n",
    "mkdir stringtie_output\n",
    "```\n",
    "Instead of simply running a line with the file names included in the command, this time we're going to define several variables and then pass the variable names to *Stringtie*. This will give the same final result, but it can be useful to use this technique when you are writing scripts which you may want to recycle later. If you want to re-run the code with a different file, it is easier to edit a script if you define all your variables in the first few lines of code. Note that the --rf option can be used in Stringtie with reads are strand-specific.\n",
    "```\n",
    "REFERENCE_ANNOTATION=\"s_cerevisiae_chrX.gff\"\n",
    "BAM_FILE=\"s_cerevisiae_chrX_sorted.bam\"\n",
    "OUTPUT_ALL=\"./stringtie_output/stringtie_all_output.gtf\"\n",
    "OUTPUT_COV=\"./stringtie_output/stringtie_cov_output.gtf\"\n",
    "\n",
    "stringtie $BAM_FILE -G $REFERENCE_ANNOTATION --rf -o $OUTPUT_ALL -A $OUTPUT_COV\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.1 Look at the output of *Stringtie*\n",
    "*Stringtie* will create some new files in the output directory you specified above. \n",
    "1. What are these files called? How many <span title=\"You may be able to use some of the commands we learned earlier, but remember: lines that are commented out (lines which begin with '#') will be counted as well if you simply use 'wc -l'...\"><b>lines</b></span> does each file contain?\n",
    "2. What is the <span title=\"Do these output files have the same columns as a SAM file, a GTF, or is it a totally different format?\"><b>format</b></span> of these files? \n",
    "\n",
    "If you want a challenge, see if you can figure out how many <span title=\"'Novel transcripts' are transcripts which were not previously known ie. not contained in the reference annotations.\"><b>novel</b></span> transcripts were assembled!"
   ]
  }
 ],
 "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}