How to determine rRNA contamination.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How to determine ribosomal RNA contamination of your libary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**_Why might you be interested in determining the amount of rRNA contamination of your library?_** Usually if you're only interested in studying polyadenylated mRNA, you would either deplete your starting material of ribosomal RNA or select for polyadenylated mRNA with oligo d(T) beads. Then, you might want to check for rRNA contamination to assess how well your depletion or selection worked.\n",
    "\n",
    "Another possibility could be that you have very limited RNA and therefore choose to avoid ribodepletion/polyA selection because these methods require a lot more starting material than you have. Then you might want to determine the amount of rRNA contamination to see if you can get away with this :)\n",
    "\n",
    "In this tutorial we will explore how to determine rRNA contamination of a library. There are multiple ways to do this and we will focus on one of them (the simpliest one given what we've already done in Module 1 and 2 in-class exercises) and touch on a few others (that may better apply to your data processing pipeline)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Method 1: Count how many reads mapped to rRNA regions with featureCounts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In Module 1 we used featureCounts to count the number of reads that mapped to each gene. Given that we already produced and [trimmed](https://github.com/biom262/biom262-2019/blob/master/Module_2/Notebooks/06_Analysis_workflows.ipynb) our counts matrix, the easiest way to determine rRNA contamination of our libraries is by manipulating this trimmed counts matrix. First let's remind ourselves of what our count matrix looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th scope=col>Geneid</th><th scope=col>Length</th><th scope=col>X.home.ucsd.train01.biom262_2019.Module_1.all_bams.NT_shRNA_hepg2_rep1_Aligned.sorted.out.bam</th><th scope=col>X.home.ucsd.train01.biom262_2019.Module_1.all_bams.NT_shRNA_hepg2_rep2_Aligned.sorted.out.bam</th><th scope=col>X.home.ucsd.train01.biom262_2019.Module_1.all_bams.TARDBP_shRNA_hepg2_rep1_Aligned.sorted.out.bam</th><th scope=col>X.home.ucsd.train01.biom262_2019.Module_1.all_bams.TARDBP_shRNA_hepg2_rep2_Aligned.sorted.out.bam</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><td>ENSG00000223972.5_2</td><td>1735               </td><td> 0                 </td><td> 0                 </td><td> 0                 </td><td> 0                 </td></tr>\n",
       "\t<tr><td>ENSG00000227232.5_2</td><td>1351               </td><td>27                 </td><td>42                 </td><td>40                 </td><td>35                 </td></tr>\n",
       "\t<tr><td>ENSG00000243485.5_5</td><td>1021               </td><td> 0                 </td><td> 0                 </td><td> 0                 </td><td> 0                 </td></tr>\n",
       "\t<tr><td>ENSG00000237613.2_2</td><td>1219               </td><td> 0                 </td><td> 0                 </td><td> 0                 </td><td> 0                 </td></tr>\n",
       "\t<tr><td>ENSG00000268020.3_4</td><td> 840               </td><td> 0                 </td><td> 0                 </td><td> 0                 </td><td> 0                 </td></tr>\n",
       "\t<tr><td>ENSG00000240361.2_3</td><td>1414               </td><td> 0                 </td><td> 0                 </td><td> 0                 </td><td> 0                 </td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|llllll}\n",
       " Geneid & Length & X.home.ucsd.train01.biom262\\_2019.Module\\_1.all\\_bams.NT\\_shRNA\\_hepg2\\_rep1\\_Aligned.sorted.out.bam & X.home.ucsd.train01.biom262\\_2019.Module\\_1.all\\_bams.NT\\_shRNA\\_hepg2\\_rep2\\_Aligned.sorted.out.bam & X.home.ucsd.train01.biom262\\_2019.Module\\_1.all\\_bams.TARDBP\\_shRNA\\_hepg2\\_rep1\\_Aligned.sorted.out.bam & X.home.ucsd.train01.biom262\\_2019.Module\\_1.all\\_bams.TARDBP\\_shRNA\\_hepg2\\_rep2\\_Aligned.sorted.out.bam\\\\\n",
       "\\hline\n",
       "\t ENSG00000223972.5\\_2 & 1735                  &  0                    &  0                    &  0                    &  0                   \\\\\n",
       "\t ENSG00000227232.5\\_2 & 1351                  & 27                    & 42                    & 40                    & 35                   \\\\\n",
       "\t ENSG00000243485.5\\_5 & 1021                  &  0                    &  0                    &  0                    &  0                   \\\\\n",
       "\t ENSG00000237613.2\\_2 & 1219                  &  0                    &  0                    &  0                    &  0                   \\\\\n",
       "\t ENSG00000268020.3\\_4 &  840                  &  0                    &  0                    &  0                    &  0                   \\\\\n",
       "\t ENSG00000240361.2\\_3 & 1414                  &  0                    &  0                    &  0                    &  0                   \\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "Geneid | Length | X.home.ucsd.train01.biom262_2019.Module_1.all_bams.NT_shRNA_hepg2_rep1_Aligned.sorted.out.bam | X.home.ucsd.train01.biom262_2019.Module_1.all_bams.NT_shRNA_hepg2_rep2_Aligned.sorted.out.bam | X.home.ucsd.train01.biom262_2019.Module_1.all_bams.TARDBP_shRNA_hepg2_rep1_Aligned.sorted.out.bam | X.home.ucsd.train01.biom262_2019.Module_1.all_bams.TARDBP_shRNA_hepg2_rep2_Aligned.sorted.out.bam | \n",
       "|---|---|---|---|---|---|\n",
       "| ENSG00000223972.5_2 | 1735                |  0                  |  0                  |  0                  |  0                  | \n",
       "| ENSG00000227232.5_2 | 1351                | 27                  | 42                  | 40                  | 35                  | \n",
       "| ENSG00000243485.5_5 | 1021                |  0                  |  0                  |  0                  |  0                  | \n",
       "| ENSG00000237613.2_2 | 1219                |  0                  |  0                  |  0                  |  0                  | \n",
       "| ENSG00000268020.3_4 |  840                |  0                  |  0                  |  0                  |  0                  | \n",
       "| ENSG00000240361.2_3 | 1414                |  0                  |  0                  |  0                  |  0                  | \n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "  Geneid              Length\n",
       "1 ENSG00000223972.5_2 1735  \n",
       "2 ENSG00000227232.5_2 1351  \n",
       "3 ENSG00000243485.5_5 1021  \n",
       "4 ENSG00000237613.2_2 1219  \n",
       "5 ENSG00000268020.3_4  840  \n",
       "6 ENSG00000240361.2_3 1414  \n",
       "  X.home.ucsd.train01.biom262_2019.Module_1.all_bams.NT_shRNA_hepg2_rep1_Aligned.sorted.out.bam\n",
       "1  0                                                                                           \n",
       "2 27                                                                                           \n",
       "3  0                                                                                           \n",
       "4  0                                                                                           \n",
       "5  0                                                                                           \n",
       "6  0                                                                                           \n",
       "  X.home.ucsd.train01.biom262_2019.Module_1.all_bams.NT_shRNA_hepg2_rep2_Aligned.sorted.out.bam\n",
       "1  0                                                                                           \n",
       "2 42                                                                                           \n",
       "3  0                                                                                           \n",
       "4  0                                                                                           \n",
       "5  0                                                                                           \n",
       "6  0                                                                                           \n",
       "  X.home.ucsd.train01.biom262_2019.Module_1.all_bams.TARDBP_shRNA_hepg2_rep1_Aligned.sorted.out.bam\n",
       "1  0                                                                                               \n",
       "2 40                                                                                               \n",
       "3  0                                                                                               \n",
       "4  0                                                                                               \n",
       "5  0                                                                                               \n",
       "6  0                                                                                               \n",
       "  X.home.ucsd.train01.biom262_2019.Module_1.all_bams.TARDBP_shRNA_hepg2_rep2_Aligned.sorted.out.bam\n",
       "1  0                                                                                               \n",
       "2 35                                                                                               \n",
       "3  0                                                                                               \n",
       "4  0                                                                                               \n",
       "5  0                                                                                               \n",
       "6  0                                                                                               "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# In Rstudio:\n",
    "# Read in trimmed counts matrix and take a look at it\n",
    "counts.matrix <- read.csv(\"~/biom262_2019/Module_1/featurecounts/featureCounts_trimmed.txt\", sep = \"\\t\")\n",
    "head(counts.matrix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For each gene ID, we have counts of how many reads mapped to it for each ouf our 4 samples. To determine rRNA contamination we will look at the portion of reads that mapped to rRNA regions. First we'll create a list of gene IDs for rRNA using the [Gencode annotation file](https://www.gencodegenes.org/human/release_29lift37.html) that we already downloaded in Module 1. If you take a look at the format of [gtf annotation files](https://uswest.ensembl.org/info/website/upload/gff.html), you'll notice that there are 9 tab separated fields and the 9th field contains semicolon separated \"attribute\" information including gene ID and gene type. Let's trim this annotation file down to just the gene IDs of rRNA gene type."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# On the command line:\n",
    "# Cut out the 9th \"attribute\" field\n",
    "tail -n +5 ~/scratch/annotations/hg19/gencode.v29lift37.annotation.gtf | cut -f9 > ~/gencode.v29lift37.annotation.attribute.gtf\n",
    "# Note that the first 4 lines of the annotation file gives some superfluous information so we ignore those."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The Gencode annotation file is almost 1.5 Gb in size, so we need an efficient way to extract \n",
    "# gene_ids associated with ribosomal RNA without having to load the full file into memory of python or R for example.\n",
    "# The awk utility provides for an easy way to do trim the large file down to just gene IDs of rRNA gene type\n",
    "awk -v RS=';' -v OFS=',' -F' ' '{if ($1==\"gene_id\") gi=$2; if($1==\"gene_type\") gt=$2}{if (gi && gt && gt ~ '/rRNA/') print gi,gt}' ~/gencode.v29lift37.annotation.attribute.gtf > rrna.geneIDs.txt\n",
    "\n",
    "# awk expression explanation:\n",
    "# -v RS=';'   # split the key-value pair records in the gtf file by semicolon\n",
    "# -v OFS=','  # use comma as the separator for the output\n",
    "# -F' '       # set space as the field separator between keys and values\n",
    "# '{if ($1==\"gene_id\") gi=$2}'   # if the left side of the match (the key) is equal to gene_id, then set the right side of the match (value) to a new variable called gi\n",
    "# 'if($1==\"gene_type\") gt=$2}'   # do the same matching for gene_type\n",
    "# {if (gi && gt && gt ~ '/rRNA/') print gi,gt}   # if both gene_id and gene_type are found in a record and gene_type matched to rRNA print both"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have a list of all the gene IDs for rRNA, we use R to easily calculate the percentage of rRNA contamination for each of our samples. I used the dplyr and reshape2 packages which you'll need to install if you haven't already."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Install dplyr and reshape2 if you haven't already\n",
    "install.packages(\"dplyr\")\n",
    "install.packages(\"reshape2\")\n",
    "\n",
    "# Read in rRNA gene IDs\n",
    "rRNA.IDs <- read.csv(\"~/rrna.geneIDs.txt\", header = FALSE)\n",
    "\n",
    "# Create a flag column to distinguish contaminant rRNA gene IDs\n",
    "counts.matrix$is_contam = counts.matrix$Geneid %in% rRNA.IDs$V1\n",
    "\n",
    "# Make sample a variable column\n",
    "count.mat.molt = melt(counts.matrix, id.vars = c('Geneid', 'Length', 'is_contam'))\n",
    "\n",
    "# Summarize percent rRNA contamination for each sample\n",
    "contam_smry = ddply(count.mat.molt, .(variable), summarise, \n",
    "                    per_contam = sum(value[is_contam])*100/sum(value))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The rRNA contamination of our samples is quite minimal (~ 4.3% for each sample). This is great! And expected because if we look closely at the [data we downloaded from Encode](https://www.encodeproject.org/experiments/ENCSR527QNC/), these libaries were prepared from polyadenylated mRNA enriched starting material (although the specific method is not noted)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Method 1a: featureCounts with rRNA only annotation file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall that we downloaded the comprehensive Gencode annotation file. What if our annotation file was not comprehensive and did not include ribosomal RNA? We could either get the comprehensive annotation file, use that to create another counts matrix with featureCounts, and do what we just did above, or we could use featureCounts with an annotation file containing only rRNA. Is there any advantage to one method or the other? Not really, but the latter might be less computationally intensive. Let's quickly go over how we would do this.\n",
    "\n",
    "First we need to make an rRNA.gtf annotation file. [Here](http://onetipperday.sterding.com/2012/02/get-rrnagtf-file-from-ucsc-table.html) are detailed step-by-step instructions for how to do this. Note that you can also include other small RNAs (e.g. tRNA, snRNA) if you'd like.\n",
    "\n",
    "Then we would use featureCounts to count how many reads aligned with features in this annotation file (just modify the script you wrote for the Module 1 homework). If we quickly look at the .err file, the percent \"Successfully assigned fragments\" should tell us the percent rRNA contamination."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Method 2: bedtools intersect"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another way to determine the amount of rRNA contamination is to use Bedtools intersect. Please see [this notebook](https://gist.github.com/byee4/4d29f1dae08581aed1e8f5bec71218e2) for detailed example code. Although this example determines overlap between CLIP peaks and CDS and 3'UTR regions, it can be modified to determine overlap between our reads and rRNA regions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Method 3: Filter rRNA contaminant reads by mapping to a \"contaminant genome\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**_What if we wanted to not only determine rRNA contamination but also filter out contaminating reads?_** We could do this by first mapping our reads to only rRNA regions and then filtering out the reads that mapped. And to figure out the percent rRNA contamination, we could just compare the number of reads that mapped to the contaminant genome with the original number of reads. Note that you could modify this approach to filter out any contaminants you'd like.\n",
    "\n",
    "1.) Specifically, we would first generate STAR genome index with only rRNA fasta files. Available [here](ftp://ftp.ebi.ac.uk/pub/databases/Rfam/CURRENT/fasta_files/) with the following Rfam accession IDs:  \n",
    "\n",
    "**5S ribosomal RNA**: RF00001  \n",
    "**5.8S ribosomal RNA**: RF00002  \n",
    "**Eukaryotic small subunit ribosomal RNA**: RF01960  \n",
    "**Eukaryotic large subunit ribosomal RNA**: RF02543  \n",
    "\n",
    "2.) Then we would map our reads to this \"rRNA contaminant genome\" (see Module 1 Notebook 4 for example scripts for steps 1 and 2). Any reads that map to it would be rRNA contaminants.\n",
    "\n",
    "3.) Finally we would filter out these contaminant reads. We could do this by creating a list of IDs of the reads that mapped to this contaminant genome and using [seqtk](https://github.com/lh3/seqtk) subseq to filter them from our fastq files."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Method 4: Filter rRNA repeat with BBDuk"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, perhaps a simpler way to filter out rRNA repeats is with [BBDuk](https://jgi.doe.gov/data-and-tools/bbtools/bb-tools-user-guide/bbduk-guide/), which will filter out repeat sequences or k-mers. For more details on how to do this see [here](https://www.biostars.org/p/159959/#159960), an example script is given below, and [here](https://www.ncbi.nlm.nih.gov/nuccore/555853?report=fasta) is the human ribosomal DNA complete repeating unit. The .err file should give you the number and percent contaminant reads."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bbduk.sh in=[YOURPATH]fastq.gz outm=ribo.fa outu=nonribo.fa k=31 ref=rRNA.fasta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see there are a lot of different ways to do roughly the same thing! Hope this tutorial was helpful!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}