Counting Transcripts - Pipeline B.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Counting Transcripts - Pipeline B\n",
    "\n",
    "> This is part three of the Counting Transcripts series. We will cover another popular pipeline: Trimmomatic/killisto. We will continue to use the data downloaded in [part 1]() of the series. This notebook should be located in the same folder as the previous notebooks.\n",
    "\n",
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "###IMPORTANT: if you ever restart this notebook, you MUST rerun this code cell\n",
    "import os\n",
    "\n",
    "#feel free to change this location, because we will be working in this folder for the entire tutorial\n",
    "os.environ['WORKDIR'] = './data' "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pipeline B: Trimmomatic/killisto\n",
    "**2a.**This pipeline will also use [Trimmomatic](http://www.usadellab.org/cms/?page=trimmomatic) to clean the reads.\n",
    "If you have walked through part 2 of the series, you can skip this step.\n",
    "If you still need to clean the reads, see instructions for using [Trimmomatic](http://www.usadellab.org/cms/?page=trimmomatic) in [pipeline A](), then return here.\n",
    "\n",
    "[Kallisto](https://pachterlab.github.io/kallisto/) is a program that appeared in [2015](https://liorpachter.wordpress.com/2015/05/10/near-optimal-rna-seq-quantification-with-kallisto/) developed in Lior Patcher's group. It is based on _pseudo-alignment_ using [k-mers](https://en.wikipedia.org/wiki/K-mer) instead of reads. It is fast, accurate, and robust. And simple to use. All you need it to create is an index of your transcriptome (or any target reference sequence) and you will be ready to apply it on your FASTQ files. This [video](https://www.youtube.com/watch?v=94wphB3GKBM) has a great explanation of how kallisto differs from other mapping tools.\n",
    "\n",
    "**2b.** The first step is to create a kallisto index, and we can use the same transcriptome FASTA file from pipeline A."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!kallisto index -i $WORKDIR/kallisto/transcripts.idx $WORKDIR/microbesOnline/4932.transcriptomes.fasta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**2c.** Kallisto combines the mapping and quantifying steps into one, so unlike pipeline A which first uses STAR to map the reads to a reference genome sequence and then uses cufflinks to quantify the aligned reads, we can run just one command to get a file representing levels of gene expression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!kallisto quant -i $WORKDIR/kallisto/transcripts.idx -o $WORKDIR/kallisto/output \\\n",
    "    $WORKDIR/cleanFASTQ/SRR5511057.forward_paired.fastq $WORKDIR/cleanFASTQ/SRR5511057.reverse_paired.fastq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once the above command has finished, you should see a new folder called `output`. If you look inside, there should be a file called `abundance.tsv` with five columns:\n",
    " 1. The first column is `target_id`, which is the identifier of the transcript from the FASTA transcriptome we downloaded. \n",
    " 2. The second column is `length`, which is the length of the transcript in base pairs. \n",
    " 3. The third column is `eff_length`, which stands for effective length. This is the length of the original transcript that matches with reads. \n",
    " 4. The fourth column is `est_counts`, which stands for estimated counts. This is almost the same as _counts_, but accounts for the amount of bias in the experiment. \n",
    " 5. This fifth column is `tpm`, which is the metric for the level of gene expression.\n",
    " \n",
    "For more information about these expression units, visit [this blog post](https://haroldpimentel.wordpress.com/2014/05/08/what-the-fpkm-a-review-rna-seq-expression-units/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Check Your Understanding\n",
    "\n",
    "1. What would we expect to be the most expressed transcript in yeast to be (hint: what do all cells need)?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Answers\n",
    "\n",
    "1. Since we know that ribosomes, the crucial component of the cell that synthesizes proteins, are made out of RNA and proteins, we would expect that a ribosomal subunit to be expressed the most. If you sort the transcripts by tpm value, we find that the highest expressed transcript is VIMSS6581441 with a TPM value of 21603.60. According to microbesOnline, VIMSS6581441 is \"RPL41B Ribosomal protein L47 of the large (60S) ribosomal subunit, identical to Rpl41Ap and has similarity to rat L41 ribosomal protein; comprised of only 25 amino acids; rpl41a rpl41b double null mutant is viable\"."
   ]
  }
 ],
 "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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}