RNAseq analysis/RNAseq_tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# RNA-seq tutorial\n",
    "*by Roseric Azondekon*\n",
    "\n",
    "**03/26/2019**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this tutorial, we show you how to download raw sequence data from the European instance of the SRA, which can\n",
    "be accessed via https://www.ebi.ac.uk/ena. At ENA, the sequencing reads are directly available in FASTQ or SRA\n",
    "formats, which will be explained below.\n",
    "\n",
    "For this tutorial, we need `FastQC`, `multiQC`, the `SRA toolkit`, the `subread` <a href=\"https://sourceforge.net/projects/subread/files/subread-1.6.4/\" target=\"_blank\">package</a>, a powerful suite of tools designed to interact with SAM and BAM files called `samtools`, `salmon`, and `STAR` installed and referenced in the environment variable `PATH`. Let's first check if this requirement is met:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "fastqc --version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "multiqc --version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fastq-dump --version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "samtools --version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "salmon -v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "STAR --version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Checking if the package is properly installed\n",
    "featureCounts -v"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If at least one of the above commands produces an error, please, check your installation of the tool and try again.\n",
    "\n",
    "Now let's create a working directory for our RNA-seq project."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mkdir -p tuto && cd tuto"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Data Download"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To download a set of SRA files:\n",
    "1. Go to https://www.ebi.ac.uk/ena.\n",
    "2. Search for the accession number of the project, e.g., SRP053046 (should be indicated in the published\n",
    "paper).\n",
    "3. There are several ways to start the download, here we show you how to do it through the command line interface on GNU/Linux.\n",
    "    - copy the link’s address of the \"SRA files\" column (right mouse click), go to the command line, move to the target directory, type: `wget < link copied from the ENA website >`\n",
    "    - If there are many samples as it is the case for the project referenced here (accession number: SRP053046), you can download the summary of the sample information from ENA by right-clicking on “TEXT\" and copying the link location."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's download the file from the link copied earlier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wget -O all_samples.txt \"https://www.ebi.ac.uk/ena/data/warehouse/filereport?\\\n",
    "accession=PRJNA274258&result=read_run&fields=study_accession,\\\n",
    "sample_accession,secondary_sample_accession,experiment_accession,\\\n",
    "run_accession,tax_id,scientific_name,instrument_model,library_layout,fastq_ftp,\\\n",
    "fastq_galaxy,submitted_ftp,submitted_galaxy,sra_ftp,sra_galaxy,cram_index_ftp,\\\n",
    "cram_index_galaxy&download=txt\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may try to open the `all_samples.txt` file with LibreOffice or Excel to view it.\n",
    "For this project, we are only interested in the paired-end first 9 RNA-seq samples. Since the first line in `all_samples.txt` contains the header, we will generate another file containing only the first 10 lines of `all_samples.txt` with the following command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sed '1d' all_samples.txt > all_samples.txt2\n",
    "head -9 all_samples.txt2 > samples.txt\n",
    "rm all_samples.txt2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's create a new folder for our SRA files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mkdir -p sra_files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "According to https://www.ncbi.nlm.nih.gov/books/NBK158899/, the FTP root to download files from NCBI is ftp://ftp-trace.ncbi.nih.gov/ and the remainder path follow the specific pattern `/sra/sra-instant/reads/ByRun/sra/{SRR|ERR|DRR}/<first 6 characters of accession>/<accession>/<accession>.sra`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the accession number for the SRA files are located in the 5th column \"Run accession\" in `all_samples.txt`. We proceed to the download of the SRA files of the samples listed in `samples.txt` with the following code: (**Attention: The download may take a long time!**)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The command below may take too long to download.\n",
    "cut -f5 samples.txt | xargs -i sh -c \\\n",
    "'v={}; FTPROOT=ftp://ftp-trace.ncbi.nih.gov/;\\\n",
    "REM=sra/sra-instant/reads/ByRun/sra/;\\\n",
    "url=${FTPROOT}${REM}${v:0:3}/${v:0:6}/${v}/${v}.sra;\\\n",
    "wget $url -P sra_files'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Converting SRA files to FASTQ files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that the download is complete, let's convert the SRA files into FASTQ files with the following command:  (**Attention: This may take a long time!**)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cut -f5 samples.txt | xargs -i sh -c \\\n",
    "'v={}; fastq-dump --outdir fastq/${v} --gzip\\\n",
    "--skip-technical --split-3 sra_files/${v}.sra'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This should be the file structure of your working directory up to this point:\n",
    "\n",
    "```shell\n",
    ".\n",
    "├── all_samples.txt\n",
    "├── samples.txt\n",
    "├── fastq\n",
    "│   ├── SRR1784137\n",
    "│   │   ├── SRR1784137_1.fastq.gz\n",
    "│   │   └── SRR1784137_2.fastq.gz\n",
    "│   ├── SRR1784138\n",
    "│   │   ├── SRR1784138_1.fastq.gz\n",
    "│   │   └── SRR1784138_2.fastq.gz\n",
    "│   ├── SRR1784139\n",
    "│   │   ├── SRR1784139_1.fastq.gz\n",
    "│   │   └── SRR1784139_2.fastq.gz\n",
    "│   ├── SRR1784140\n",
    "│   │   ├── SRR1784140_1.fastq.gz\n",
    "│   │   └── SRR1784140_2.fastq.gz\n",
    "│   ├── SRR1784141\n",
    "│   │   ├── SRR1784141_1.fastq.gz\n",
    "│   │   └── SRR1784141_2.fastq.gz\n",
    "│   ├── SRR1784142\n",
    "│   │   ├── SRR1784142_1.fastq.gz\n",
    "│   │   └── SRR1784142_2.fastq.gz\n",
    "│   ├── SRR1784143\n",
    "│   │   ├── SRR1784143_1.fastq.gz\n",
    "│   │   └── SRR1784143_2.fastq.gz\n",
    "│   ├── SRR1784144\n",
    "│   │   ├── SRR1784144_1.fastq.gz\n",
    "│   │   └── SRR1784144_2.fastq.gz\n",
    "│   └── SRR1784145\n",
    "│       ├── SRR1784145_1.fastq.gz\n",
    "│       └── SRR1784145_2.fastq.gz\n",
    "└── sra_files\n",
    "    ├── SRR1784137\n",
    "    ├── SRR1784137.sra\n",
    "    ├── SRR1784138\n",
    "    ├── SRR1784138.1\n",
    "    ├── SRR1784138.sra\n",
    "    ├── SRR1784139\n",
    "    ├── SRR1784139.sra\n",
    "    ├── SRR1784140.sra\n",
    "    ├── SRR1784141.sra\n",
    "    ├── SRR1784142.sra\n",
    "    ├── SRR1784143.sra\n",
    "    ├── SRR1784144.sra\n",
    "    └── SRR1784145.sra\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Quality Control of the FASTQ files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Up to this point, we have all our RNA-seq FASTQ files ready for Quality Control (QC) check. This is done with the `fastqc` tools developed by the <a href=\"https://www.bioinformatics.babraham.ac.uk/projects/fastqc/\" target=\"_blank\">Babraham Institute</a>. Run the following command to perform QC check for all the samples: (**This may take some time!**)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cut -f5 samples.txt | xargs -i sh -c \\\n",
    "'v={}; mkdir -p fastqc_reports/${v}; fastqc\\\n",
    "fastq/${v}/*fastq.gz -o fastqc_reports/${v}'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, let's summarize the QC reports (for all the samples) into one unique report using `multiqc`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "multiqc fastqc_reports --dirs -o multiQC_report/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's examine the summary `multiqc` report either by double-clicking on `multiQC_report/multiqc_report.html` or by executing the following code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xdg-open multiQC_report/multiqc_report.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Read Alignment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to identify the transcripts that are present in a specific sample, the genomic origin of the sequenced\n",
    "cDNA fragments must be determined. The assignment of sequencing reads to the most likely locus of origin\n",
    "is called read alignment or mapping and it is a crucial step in most types of high-throughput sequencing\n",
    "experiments.\n",
    "\n",
    "The general challenge of short read alignment is to map millions of reads accurately and in a reasonable time,\n",
    "despite the presence of sequencing errors, genomic variation and repetitive elements. The different alignment\n",
    "programs employ various strategies that are meant to speed up the process (e.g., by indexing the reference\n",
    "genome) and find a balance between mapping fidelity and error tolerance."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.1. Reference genomes and annotation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Genome sequences and annotation are often generated by consortia such as (mod)ENCODE, The Mouse\n",
    "Genome Project, The Berkeley Drosophila Genome Project, and many more. The results of these efforts can\n",
    "either be downloaded from individual websites set up by the respective consortia or from more comprehensive data bases such as the one hosted by the University of California, Santa Cruz (<a href=\"https://genome.ucsc.edu/\" target=\"_blank\">UCSC</a>) or the European genome resource (<a href=\"http://www.ensembl.org\" target=\"_blank\">Ensembl</a>).\n",
    "\n",
    "Reference sequences are usually stored in plain text FASTA files that can either be compressed with the\n",
    "generic gzip command. The annotation file is often stored as a GTF (Gene Transfer Format) despite the availability of several other file formats.\n",
    "\n",
    "Both reference sequences and GTF annotation file can be obtained either from <a href=\"https://www.ncbi.nlm.nih.gov/genome/51\" target=\"_blank\">NCBI</a>, <a href=\"https://www.ensembl.org/info/data/ftp/index.html\" target=\"_blank\">ENSEMBL</a> or <a href=\"http://hgdownload.soe.ucsc.edu/downloads.html#human\" target=\"_blank\">UCSC Genome Browser</a>.\n",
    "\n",
    "While it is usually faster to align against a cDNA reference sequence (cDNA), in this tutorial, we will align the reads against the genome (DNA) reference sequences. We obtain both the genome refernce sequences and our gene annotation files from <a href=\"https://www.ensembl.org/info/data/ftp/index.html\" target=\"_blank\">ENSEMBL</a>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Download the latest human genome\n",
    "wget -P reference \\\n",
    "ftp://ftp.ensembl.org/pub/release-93/fasta/homo_sapiens\\\n",
    "/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Decompress genome sequence file and rename it\n",
    "gzip -dk < reference/Homo_sapiens.GRCh38.dna.\\\n",
    "primary_assembly.fa.gz > reference/human_genome.fa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Download GTF annotation file\n",
    "wget -P annotation ftp://ftp.ensembl.org/pub/release-95\\\n",
    "/gtf/homo_sapiens/Homo_sapiens.GRCh38.95.gtf.gz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Decompress gene annotation file and rename it\n",
    "gzip -dk < annotation/Homo_sapiens.GRCh38.95.gtf.gz\\\n",
    "> annotation/gene_annotation.gtf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.2. Aligning reads using STAR"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.2.1. Generate genome index"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**This step has to be done only once per genome type (and alignment program)**. The index files will comprise the genome sequence, suffix arrays (i.e., tables of k-mers), chromosome names and lengths, splice junctions coordinates, and information about the genes (e.g. the strand). Therefore, the main input for this step encompasses the reference genome and an annotation file.\n",
    "\n",
    "**The index creation may take a long time!**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a directory to store the index in\n",
    "mkdir -p STARindex"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run STAR in \"genomeGenerate\" mode with 10 threads\n",
    "# Feel free to change the 'runThreadN' parameter \n",
    "# depending on the number of cores available on your computer\n",
    "STAR --runThreadN 10 \\\n",
    "--runMode genomeGenerate --genomeDir STARindex \\\n",
    "--genomeFastaFiles reference/human_genome.fa \\\n",
    "--genomeChrBinNbits 15 \\\n",
    "--sjdbGTFfile annotation/gene_annotation.gtf \\\n",
    "--sjdbOverhang 49"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.2.2. Alignment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This step has to be done for each individual FASTQ file.\n",
    "\n",
    "**This step may take a long time!**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a directory to store the alignment files\n",
    "mkdir -p alignment_STAR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# execute STAR in the runMode \"alignReads\"\n",
    "cut -f5 samples.txt | xargs -i sh -c \\\n",
    "'v={}; STAR --genomeDir STARindex \\\n",
    "--readFilesIn fastq/${v}/*fastq.gz \\\n",
    "--readFilesCommand zcat \\\n",
    "--outFileNamePrefix alignment_STAR/${v}\\\n",
    "--outFilterMultimapNmax 1 \\\n",
    "--outReadsUnmapped Fastx \\\n",
    "--outSAMtype BAM SortedByCoordinate \\\n",
    "--twopassMode Basic \\\n",
    "--runThreadN 10'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.2.3. BAM file indexing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most downstream applications like the <a href=\"https://software.broadinstitute.org/software/igv/\" target=\"_blank\">Integrative Genomics Viewer (IGV)</a> will require a .BAM.BAI file together with every BAM file to quickly access the BAM files without having to load them into memory.\n",
    "\n",
    "To obtain these index, let's now index all the BAM files within the 'alignment_STAR' folder using `samtools`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in alignment_STAR/*; do\n",
    "    if [ \"${i}\" != \"${i%.bam}\" ];then\n",
    "        samtools index ${i}\n",
    "    fi\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.3. Read alignment assessment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are numerous ways to do basic checks of the alignment success. An alignment of RNA-seq reads is\n",
    "usually considered to have succeeded if the mapping rate is >70%.\n",
    "\n",
    "The very first QC of aligned reads should be to generally check the aligner’s output. We can perform a basic assessment of the read alignment for all our samples with using the `samtools flagstat` command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in alignment_STAR/*; do\n",
    "    if [ \"${i}\" != \"${i%.bam}\" ];then\n",
    "        echo ${i:15:10}\n",
    "        samtools flagstat ${i}\n",
    "        echo '############################################################'\n",
    "    fi\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also inspect the final log file generated by STAR with the following command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in alignment_STAR/*; do\n",
    "    if [ \"${i}\" != \"${i%.final.out}\" ];then\n",
    "        echo ${i:15:10}\n",
    "        cat ${i}\n",
    "        echo\n",
    "        echo '############################################################'\n",
    "        echo \n",
    "    fi\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5. Read Quantification"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5.1. Using STAR output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To compare the expression of single genes between different conditions, an essential step is the quantification\n",
    "of reads per gene. The most popular tools for gene quantification are `htseq-count` and `featureCounts`. `featureCounts` is provided by the `subread` package.\n",
    "\n",
    "The `featureCounts` tool calls a hit if any overlap (1 bp or more) is found between the read and a feature and provides the option to either exclude multi-overlap reads or to count them for each feature that is overlapped.\n",
    "\n",
    "The nature and lengths of the reads, gene expression quantification will be strongly affected by the underlying gene annotation file that is supplied to the quantification programs.\n",
    "\n",
    "In the following command, `featureCounts` is used to count the number of reads overlapping with genes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a folder for read counts\n",
    "mkdir -p read_counts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Count features using 10 threads ('-T 10')\n",
    "featureCounts -a annotation/gene_annotation.gtf\\\n",
    "              -o read_counts/featureCounts_results.txt\\\n",
    "                 alignment_STAR/*bam -T 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5.2. Using the `salmon` pseudoaligner"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`salmon` is a probabilistic RNA-seq quantification program. It pseudoailigns reads to a reference sequence, producing a list of transcripts that are compatible with each read while avoiding alignment of individual bases.\n",
    "\n",
    "The program circumvents the need for large alignment files, thus reducing storage needs while increasing speed and enabling the processing of large numbers of samples on modest computational resources. Pseudoaligners like `salmon` estimate **transcript-level counts (not gene-level counts)**. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.2.1. Transcriptome Indexing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now download the Homo sapiens transcriptome reference from <a href=\"https://www.ncbi.nlm.nih.gov/genome\" target=\"_blank\">NCBI</a>. Transcriptome reference sequence file can also be downloaded from <a href=\"https://www.ensembl.org/info/data/ftp/index.html\" target=\"_blank\">ENSEMBL</a>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wget -P reference ftp://ftp.ncbi.nlm.nih.gov/refseq/H_sapiens/annotation\\\n",
    "/GRCh38_latest/refseq_identifiers/GRCh38_latest_rna.fna.gz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Decompress transcriptome sequence file and rename it\n",
    "gzip -dk < reference/GRCh38_latest_rna.fna.gz > reference/human_transcriptome.fna"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, let's build an index on our transcriptome (**this may take some time!**)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if [ ! -d salmon_index ] ; then\n",
    "    salmon index -t reference/human_transcriptome.fna -i salmon_index\n",
    "fi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.2.2. Quantifying reads using `salmon`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have our index built, we are ready to quantify our samples using the following command (**this may take some time!**):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read quantification using 8 threads (-p 8)\n",
    "cut -f5 samples.txt | xargs -i sh -c \\\n",
    "'v={}; salmon quant -i salmon_index -l A \\\n",
    "         -1 fastq/${v}/${v}_1.fastq.gz \\\n",
    "         -2 fastq/${v}/${v}_2.fastq.gz \\\n",
    "         -p 8 -o quants/${v}_quant'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The read counts file `featureCounts_results.txt` can be imported in `R` for downstream data analysis. Similarly, we can import transcript-level estimates into `R` from the quantification files generated by `salmon`.\n",
    "\n",
    "In another tutorial, we will show how to process `featureCounts_results.txt` and the transcript level files from `salmon` for Differential Gene Expression Analysis using `DESeq2`, `edgeR`, or `limma-voom` in `R`."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Bash",
   "language": "bash",
   "name": "bash"
  },
  "language_info": {
   "file_extension": ".sh",
   "help_links": [
    {
     "text": "MetaKernel Magics",
     "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
    }
   ],
   "mimetype": "text/x-bash",
   "name": "bash"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}