Counting Transcripts - Introduction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Counting Transcripts - An Introduction\n",
    "\n",
    "> Cell phenotype is controlled by gene expression regulation, a phenomenon discovered by Barbara McClintock in early 1950s, later described by  Francois Jacob and Jacques Monod in bacteria. Since the discovery of the entire human DNA sequence, scientists have now the knowledge and technology to quantify precisely gene expression at whole genome scale. A particularly interesting layer of regulation occurs at the transcriptome level. The transcriptome is the set of all RNA molecules in the cell. Cellular phenotype is particularly impacted by the pool of messenger RNAs (mRNAs), which results from the transcription of specific genes and serves as template for protein translation. Through the use of transcriptomics techniques, we can determine the abundance of each mRNA in a cell, therefore determining the levels of gene expression. Quantifying gene expression levels in our cells helps us understand important biological phenomena like cell differentiation and mechanisms of human disease.\n",
    "\n",
    "> There are two main transcriptomics technologies use to measure gene expression: (1) [Microarrays](https://en.wikipedia.org/wiki/DNA_microarray) is an older technique used to measure relative abundance of gene expression. We will not cover microarrays here, but this [video](https://www.youtube.com/watch?v=6ZzFihESjp0) does a great job explaining how microarrays works, and (2) RNA sequencing (RNA-seq). In this tutorial series, we will focus  on RNA-seq data analysis.\n",
    "\n",
    "\n",
    "##  RNA-seq\n",
    "\n",
    "> [RNA-seq](https://en.wikipedia.org/wiki/RNA-Seq) is a technology to quantify gene expression by determining the abundance of specific RNA sequences. There are three main steps to RNA-seq: library preparation, sequencing, and data analysis. \n",
    "\n",
    "> Let's start with the steps happening in the lab. First, we isolate mRNA. Then, using the enzyme reverse transcriptase, mRNA is converted into its complementary DNA strand, cDNA. The cDNA is then cut up into small fragments, usually a few hundred base-pairs long. [Next generation sequencing technologies](https://en.wikipedia.org/wiki/Illumina_dye_sequencing) use fluorescence-tagged nucleotides that bind to short cDNA sequences to recapitulate the sequence of the sample by sequencing, in a massively parallel manner. Each sequenced read (call *reads*) are stored into a file of [FASTQ format](https://en.wikipedia.org/wiki/FASTQ_format).\n",
    "\n",
    "> Later, as part of the data analysis pipeline, these files containing the sequence of the reads are cleaned using quality control pipelines, then mapped from their sequence to a specific genome / transcriptome / set of features.\n",
    "\n",
    "## File Structure\n",
    "\n",
    "> Organization is a very important skill for any bioinformatics project that is oftentimes overlooked. There is no standard file hierarchy to use, but we will use a file structure that is fairly common among many research scientists:\n",
    "\n",
    "<html>\n",
    "<div style=\"width: 50%;margin: 0 auto; background-color: rgb(242, 242, 242);\">\n",
    "    <code><i>project_name</i>/\n",
    "    data/\n",
    "        BAM/\n",
    "        cleanFASTQ/\n",
    "        cufflinks/\n",
    "        FASTQ/\n",
    "        kallisto/\n",
    "        microbesOnline/\n",
    "        seqData/\n",
    "        transcriptomeIndex/</code>\n",
    "</div>\n",
    "</html>\n",
    "\n",
    "> Also, we will store all downloaded software in a  folder called *software* in the home directory, ~/software, and add binaries of downloaded software to the [PATH](https://en.wikipedia.org/wiki/PATH_(variable)). Instructions for adding to your PATH on Mac can be found [here](https://coolestguidesontheplanet.com/add-shell-path-osx/).\n",
    "\n",
    "## Setup\n",
    "\n",
    "> Below, we will assign the path to the working directory as a global variable that will be used throughout the rest of the notebook."
   ]
  },
  {
   "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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#create the directory if doesn't exist already\n",
    "!mkdir -p $WORKDIR"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# RNA-Seq Data Analysis\n",
    "> In this tutorial, we will use a subset of samples from a published work, where they investigated how microbes adapt to new environments ([PMID 28854640](https://www.ncbi.nlm.nih.gov/pubmed/28854640)). Following standard practices in science, the authors made all data required to reproduce their findings, including the reads for several transcriptome samples, publically available (SRA Project [SRP106486](https://trace.ncbi.nlm.nih.gov/Traces/sra/?study=SRP106486)).\n",
    "\n",
    "First, we need to download the raw data, which in this case is a file of 500 MB. For a connection speed of around 6 MB/s, the file will download in about 2 minutes, but download time will vary depending on your internet connection. Afterwards, move the file into the $WORKDIR/seqData/ folder. \n",
    "\n",
    "Note: Access to this file may have changed. In this case, tryhttps://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos2/sra-pub-run-7/SRR5511057/SRR5511057.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!wget ftp://ftp-trace.ncbi.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR551/SRR5511057/SRR5511057.sra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "If the above command returned `wget: command not found`, means you don't have the command `wget` in your machine. Then, try installing Homebrew, run `brew install wget` in terminal, and then try the above command again."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preparing the file\n",
    "You may notice that we downloaded a `.sra` file. This file format contains the raw sequencing data and alignment information in a compressed form. We will convert this into [FASTQ](https://en.wikipedia.org/wiki/FASTQ_format) format, which is the standard file format for storing RNA-seq raw data. FASTQ is a text-based format for storing the biological sequence and its corresponding quality scores. Each read is encoded with 4 lines in strict order:\n",
    "\n",
    "* Line 1: beings with '@' character with an identifier for the sequence and an optional description.\n",
    "* Line 2: raw nucleotide sequence.\n",
    "* Line 3: begins with '+' character with the same sequence identifier and additional optional description.\n",
    "* Line 4: quality of each nucleotide.\n",
    "\n",
    "To help us convert this file into FASTQ, you will need to download [SRA Toolkit](https://www.ncbi.nlm.nih.gov/books/NBK158900/).\n",
    "Although this toolkit contains many helpful tools to help us work with sequencing data, there is only one command we need to learn for now:\n",
    "\n",
    "`fastq-dump --split-files -O <outputDir> <sraFile>`\n",
    "\n",
    "This command converts the .sra file into FASTQ format. The --split-files tag specifies whether our data was obtained using single-end sequencing or paired-end sequencing. Single-end sequencing means that each hybridized fragment is read once, while paired-end sequencing means that each hybridized fragment is read twice (once from each end). Paired-end sequencing results in twice as many reads compared to single-end sequencing, so is therefore more costly, but improves accuracy and the quality of the sequencing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!fastq-dump --split-files -O $WORKDIR/FASTQ $WORKDIR/seqData/SRR5511057.sra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our data was sequenced using paired-end sequencing, so we used the `--split-files` flag. You should see two new files in your $WORKDIR/FASTQ folder: `SRR5511057_1.fastq` and `SRR5511057_2.fastq`. If you inspect the SRR5511057_1.fastq file with the `head` command, we can see what our reads look like in FASTQ format. Here is the first read, corresponding to the first 4 lines of the file:\n",
    "\n",
    "<div style=\"border:solid 1px; padding:5px; display: inline-block\">\n",
    "@SRR5511057.1 1 length=76 <br>\n",
    "GTCCACACAAGCCAAATGGTGAAATGCAGTTCTGGGTGTGCCTAGGTATAAGGATTTGTCCTTCTTGCTGACAGCT <br>\n",
    "+SRR5511057.1 1 length=76 <br>\n",
    "AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE/EEEEEEEEEEEEEAEEEEE\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Converting Raw Data into Usable Data\n",
    "\n",
    "There are three overall steps in RNA-seq transcriptome qualification:\n",
    "1. **Cleaning reads.** Reads typically have a length of 75-250 bp (base-pairs). Due to molecular noise, not all positions of the read have the same likelihood of that event being true. For this reason, we have to clean and trim the data.\n",
    "2. **Mapping reads.** We need to associate reads to our reference features, e.g., transcripts. For this task, we need the FASTQ files and a reference annotation: (1) a set of files that define the genome, typically FASTA files and (2) a file defining where transcripts are located in the genome, typically a GFF3 file. The output of the mapping algorithms is a BAM file, which is a compressed binary file that contains the aligned reads to the reference.\n",
    "3. **Quantifying reads.** Next, we need to associate an expression value to each feature, i.e., quantify the level of expression of our transcripts from the reads alignment.\n",
    "\n",
    "There are 3 major units for quantifying reads: \n",
    "- FPKM: is useful for visualization across genes and samples because it accounts for the variability of transcript length and experimental sequencing depth. This is useful, since longer transcripts will produce more transcript fragments by chance, and in turn will have more reads associated with it since the longer the transcript.\n",
    "- TPM: is similar to FPKM, but it normalizes each gene expression to a constant derived from the sum of expression across all genes. This is helpful for comparison among samples where effective transcript lenght changes. Here's a very helpful [blog post and video](https://www.rna-seqblog.com/rpkm-fpkm-and-tpm-clearly-explained/) understanding the difference between TPM and FPKM.\n",
    "- _counts_: is the raw number of reads that match to the transcript. It is useful for statistical purposes. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Check Your Understanding\n",
    "\n",
    "1. Why are the second and fourth lines of the FASTA format the same length?\n",
    "2. Why is paired end sequencing preferred over single end sequencing?\n",
    "3. Why does fastq-dump create two files?\n",
    "4. If the sequence is 150bp apart, yet on the reference the reads are 100bp apart, does this indicate an insertion or a deletion in the genome?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Answers\n",
    "\n",
    "1. Each letter in the fourth line represents the quality of the corresponding nucleotide, so the number of quality scores must match the total number of nucleotides in the sequence.\n",
    "2. Paired end sequencing improves the accuracy of sequencing. By knowing the approximate length of the sequence and the nucleotide sequence on both ends, we can more accurately map the sequence to its location on the reference genome.\n",
    "3. Each file for each read. Reads have the same order in both files, they are paired. This [link](https://en.wikipedia.org/wiki/FASTQ_format) explains more details about the format.\n",
    "4. Insertion, for the 100bp to become 150bp, there needs to have been an insertion of ~50 bp between the two reads."
   ]
  }
 ],
 "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
}