docs/ipynb/myvariant_R_miller.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MyVariant.info and MyGene.info Use Case"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following R script demonstrates the utility of the **MyVariant.info** and **MyGene.info** R clients to annotate variants and prioritize candidate genes in patients with rare Mendelian diseases. This specific study uses data obtained from the database of phenotype and genotype (dbGaP) study.  FASTQ files generated by Ng et al for the [Miller syndrome study](http://www.ncbi.nlm.nih.gov/pubmed/19915526) were processed according to the Broad Institute’s best practices. Individual samples were aligned to the hg19 reference genome using BWA-MEM 0.7.10. Variants were called using GATK 3.3-0 HaplotypeCaller and quality scores were recalibrated using GATK VariantRecalibrator."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Initial Library Imports and Data Loading"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n",
      "Attaching package: ‘BiocGenerics’\n",
      "\n",
      "The following objects are masked from ‘package:parallel’:\n",
      "\n",
      "    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,\n",
      "    clusterExport, clusterMap, parApply, parCapply, parLapply,\n",
      "    parLapplyLB, parRapply, parSapply, parSapplyLB\n",
      "\n",
      "The following object is masked from ‘package:stats’:\n",
      "\n",
      "    xtabs\n",
      "\n",
      "The following objects are masked from ‘package:base’:\n",
      "\n",
      "    anyDuplicated, append, as.data.frame, as.vector, cbind, colnames,\n",
      "    do.call, duplicated, eval, evalq, Filter, Find, get, intersect,\n",
      "    is.unsorted, lapply, Map, mapply, match, mget, order, paste, pmax,\n",
      "    pmax.int, pmin, pmin.int, Position, rank, rbind, Reduce, rep.int,\n",
      "    rownames, sapply, setdiff, sort, table, tapply, union, unique,\n",
      "    unlist, unsplit\n",
      "\n",
      "Creating a generic function for ‘nchar’ from package ‘base’ in package ‘S4Vectors’\n",
      "\n",
      "Attaching package: ‘VariantAnnotation’\n",
      "\n",
      "The following object is masked from ‘package:base’:\n",
      "\n",
      "    tabulate\n",
      "\n",
      "Welcome to Bioconductor\n",
      "\n",
      "    Vignettes contain introductory material; view with\n",
      "    'browseVignettes()'. To cite Bioconductor, see\n",
      "    'citation(\"Biobase\")', and for packages 'citation(\"pkgname\")'.\n",
      "\n",
      "Loading required package: DBI\n",
      "\n",
      "\n",
      "Attaching package: ‘plyr’\n",
      "\n",
      "The following object is masked from ‘package:XVector’:\n",
      "\n",
      "    compact\n",
      "\n",
      "The following object is masked from ‘package:IRanges’:\n",
      "\n",
      "    desc\n",
      "\n",
      "The following object is masked from ‘package:S4Vectors’:\n",
      "\n",
      "    rename\n",
      "\n"
     ]
    }
   ],
   "source": [
    "library(myvariant, quietly=TRUE)\n",
    "library(mygene, quietly=TRUE)\n",
    "library(VariantAnnotation, quietly=TRUE)\n",
    "library(GO.db, quietly=TRUE)\n",
    "source(\"https://raw.githubusercontent.com/SuLab/myvariant.info/master/docs/ipynb/mendelian.R\")\n",
    "setwd(\"~/sulab/myvariant/vcf/recal\")\n",
    "vcf.files <- paste(getwd(), list.files(getwd()), sep=\"/\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`vcf.files` contains paths to the vcf files for each of the four patients included in this analysis.  Exome sequence data from two sibs with Miller syndrome and two unrelated affected individuals used in this vignette was provided by [Ng et al. (2010) Nature Genetics (phs000244.v1.p1)](http://www.ncbi.nlm.nih.gov/pubmed/19915526).  As this is protected information, access must be requested from dbGaP [here](http://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/study.cgi?study_id=phs000244.v1.p1) in order to run this notebook.\n",
    "\n",
    "mendelian.R defines some helper functions that are used in the analysis occurring after annotation retrieval:\n",
    "\n",
    "`replaceWith0` - replaces all NAs in a data.frame with 0.\n",
    "\n",
    "`rankByCaddScore` - for prioritizing genes by deleteriousness (scaled CADD score)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Annotating variants with MyVariant.info\n",
    "The following function reads in each output VCF file using the VariantAnnotation package available from Bioconductor. Install with `biocLite(\"VariantAnnotation\")`. `formatHgvs` (from the ***myvariant*** Bioconductor package) is a function that reads the genomic location and variant information from the VCF to create HGVS IDs which serve as a primary key for each variant. The function `getVariants` makes the queries to MyVariant.info to retrieve annotations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing  /Users/cyrusafrasiabi/recal/subject01_recalibrate_SNP_vqsr.vcf ...\n",
      "found header lines for 3 ‘fixed’ fields: ALT, QUAL, FILTER \n",
      "found header lines for 24 ‘info’ fields: AC, AF, ..., VQSLOD, culprit \n",
      "found header lines for 5 ‘geno’ fields: GT, AD, DP, GQ, PL \n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Concatenating data, please be patient.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing  /Users/cyrusafrasiabi/recal/subject02_recalibrate_SNP_vqsr.vcf ...\n",
      "found header lines for 3 ‘fixed’ fields: ALT, QUAL, FILTER \n",
      "found header lines for 24 ‘info’ fields: AC, AF, ..., VQSLOD, culprit \n",
      "found header lines for 5 ‘geno’ fields: GT, AD, DP, GQ, PL \n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Concatenating data, please be patient.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing  /Users/cyrusafrasiabi/recal/subject03_recalibrate_SNP_vqsr.vcf ...\n",
      "found header lines for 3 ‘fixed’ fields: ALT, QUAL, FILTER \n",
      "found header lines for 24 ‘info’ fields: AC, AF, ..., VQSLOD, culprit \n",
      "found header lines for 5 ‘geno’ fields: GT, AD, DP, GQ, PL \n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Concatenating data, please be patient.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing  /Users/cyrusafrasiabi/recal/subject04_recalibrate_SNP_vqsr.vcf ...\n",
      "found header lines for 3 ‘fixed’ fields: ALT, QUAL, FILTER \n",
      "found header lines for 24 ‘info’ fields: AC, AF, ..., VQSLOD, culprit \n",
      "found header lines for 5 ‘geno’ fields: GT, AD, DP, GQ, PL \n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Concatenating data, please be patient.\n"
     ]
    }
   ],
   "source": [
    "getVars <- function(vcf.file){\n",
    "  cat(paste(\"Processing \", vcf.file, \"...\\n\", sep=\" \"))\n",
    "  vcf <- readVcf(vcf.file, genome=\"hg19\")\n",
    "  vcf <- vcf[isSNV(vcf)]\n",
    "  vars <- rowRanges(vcf)\n",
    "  vars <- as(vars, \"DataFrame\")\n",
    "  vars$query <- formatHgvs(vcf, \"snp\")\n",
    "  annotations <- getVariants(vars$query, fields=c(\"dbnsfp.genename\", \"dbnsfp.1000gp1.af\", \n",
    "                                                  \"exac.af\", \"cadd.consequence\", \"cadd.phred\"), verbose=FALSE)\n",
    "  annotations[c('DP', 'FS', 'QD')] <- info(vcf)[c('DP', 'FS', 'QD')]\n",
    "  annotations <- replaceWith0(annotations)\n",
    "  annotations <- subset(annotations, !(dbnsfp.genename %in% c(\"NULL\", 0)))\n",
    "  annotations\n",
    "}\n",
    "\n",
    "vars <- lapply(vcf.files, getVars)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All genes (variants with a valid `dbnsfp.genename`) that are mutated amongst all four patients are examined.  The following function counts the number of genes in `inp` that are mutated among all four patients:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "countGenes <- function(inp) {\n",
    "    ret <- subset(data.frame(table(unlist(lapply(inp, function(i) unique(i$dbnsfp.genename))))), \n",
    "    Freq == 4)\n",
    "    cat(\"Genes remaining: \", paste(nrow(ret)))\n",
    "    ret\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Initial Number of Genes Mutated in All Patients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Genes remaining:  2441"
     ]
    }
   ],
   "source": [
    "nVars <- countGenes(vars)                                                        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Filtering Steps\n",
    "The next steps apply annotation based filters by subsetting each data.frame by those variants that meet certain criteria.\n",
    "\n",
    "### 1 - Filtering for Coverage and Strand Bias\n",
    "The first filter utilizes annotations output by HaplotypeCaller from GATK. Variants are kept that have at least 8 reads, Fisher Strand bias less than 30, and quality over depth greater than 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Genes remaining:  2308"
     ]
    }
   ],
   "source": [
    "filter1 <- lapply(vars, function(i) subset(i, DP > 8 & FS < 30 & QD > 2))\n",
    "\n",
    "nFilter1 <- countGenes(filter1) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2 - Filtering for Nonsynonymous and Splice Site Variants\n",
    "Mendelian diseases are most likely to be caused by nonsynonymous mutations. The CADD database annotates the mutation type in the field \"cadd.consequence\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Genes remaining:  1917"
     ]
    }
   ],
   "source": [
    "filter2 <- lapply(filter1, function(i) subset(i, cadd.consequence %in% c(\"NON_SYNONYMOUS\", \"STOP_GAINED\", \"STOP_LOST\", \n",
    "                                           \"CANONICAL_SPLICE\", \"SPLICE_SITE\")))\n",
    "\n",
    "nFilter2 <- countGenes(filter2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3 - Filtering for Allele Frequency Annotated by ExAC\n",
    "The third filter keeps rare variants according to the ExAC data set with allele frequency < 0.01. Rare diseases are likely caused by mutations that have not been documented yet."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Genes remaining:  18"
     ]
    }
   ],
   "source": [
    "filter3 <- lapply(filter2, function(i) subset(i, exac.af < 0.01))\n",
    "\n",
    "nFilter3 <- countGenes(filter3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4 - Filtering for Allele Frequency Annotated by 1000 Genomes Project\n",
    "The fourth filter keeps rare variants according to the 1000 Genomes Project with allele frequency < 0.01."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Genes remaining:  9"
     ]
    }
   ],
   "source": [
    "filter4 <- lapply(filter3, function(i) subset(i, sapply(dbnsfp.1000gp1.af, function(j) j < 0.01 )))\n",
    "    \n",
    "top.genes <- countGenes(filter4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5 - Filtering by GO Biological Process Annotation using MyGene.info\n",
    "Since Miller Syndrome is known to be an inborn error of metabolism, this filter keeps only genes involved in metabolic processes according to their GO biological process annotation.  To accomplish this, GO biological process annotations are pulled for each remaining gene using the **MyGene.info** R client, which can be installed from Bioconductor (`biocLite(\"mygene\")`).  Here, the `queryMany` function is used, requesting the necessary annotations using the `fields` parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Finished\n"
     ]
    }
   ],
   "source": [
    "goBP <- data.frame(queryMany(top.genes$Var1, scopes=\"symbol\", species=\"human\", fields=c(\"go.BP\", \"name\", \"MIM\", \"uniprot\")))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Bioconductor package **go.DB** is used to find all genes with a GO biological process annotation that is a descendant of `GO:0008152` - the GO id for metabolic process."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Genes remaining:  5"
     ]
    }
   ],
   "source": [
    "miller.bp <- lapply(goBP$go.BP, function(i) unlist(i$id))\n",
    "bp.ancestor <- lapply(miller.bp, function(i) sapply(i, function(j) \"GO:0008152\" %in% unlist(GOBPANCESTOR[[j]])))\n",
    "candidate.genes <- top.genes$Var1[sapply(bp.ancestor, function(i) TRUE %in% i)]\n",
    "cat(\"Genes remaining: \", length(candidate.genes))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Prioritizing genes\n",
    "The remaining five genes can be prioritized according to CADD (deleteriousness) score. `rankByCaddScore` extracts the average CADD scores of the variants in each gene and ranks in descending order."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th></th><th scope=col>gene</th><th scope=col>cadd.phred</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>1</th><td>DHODH</td><td>26.81</td></tr>\n",
       "\t<tr><th scope=row>2</th><td>CTBP2</td><td>21.385</td></tr>\n",
       "\t<tr><th scope=row>3</th><td>PIK3R3</td><td>20.7</td></tr>\n",
       "\t<tr><th scope=row>4</th><td>CDC27</td><td>18.545</td></tr>\n",
       "\t<tr><th scope=row>5</th><td>CDON</td><td>10.02</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|ll}\n",
       "  & gene & cadd.phred\\\\\n",
       "\\hline\n",
       "\t1 & DHODH & 26.81\\\\\n",
       "\t2 & CTBP2 & 21.385\\\\\n",
       "\t3 & PIK3R3 & 20.7\\\\\n",
       "\t4 & CDC27 & 18.545\\\\\n",
       "\t5 & CDON & 10.02\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/plain": [
       "    gene cadd.phred\n",
       "1  DHODH     26.810\n",
       "2  CTBP2     21.385\n",
       "3 PIK3R3     20.700\n",
       "4  CDC27     18.545\n",
       "5   CDON     10.020"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ranked <- rankByCaddScore(candidate.genes, filter4)\n",
    "ranked"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "This analysis highlights the use of the **MyVariant.info** and **MyGene.info** annotation services to narrow the candidate gene list from 2441 genes to 5 - representing a significant decrease in the burden of manual biological analysis."
   ]
  }
 ],
 "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.2.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}