COO5b.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Practical exercise sheet 5:\n",
    "# Phasing of genetic sequencing data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  PURPOSE: Obtain a VCF file where all (/as many as possible of) the variant genotypes of sample NA12878, across brca1 and brca2 genes, are phased !\n",
    "We are going to test two of the methods for phasing presented during the lectures:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.Pedigree based phasing of genetic variants (using GATK PhaseByTransmission - PBT)\n",
    "\n",
    "\n",
    "You can look at PBT documentation to see more details about what it needs and what it does here:\n",
    "https://software.broadinstitute.org/gatk/gatkdocs/3.5-0/org_broadinstitute_gatk_tools_walkers_phasing_PhaseByTransmission.php\n",
    "\n",
    "<span style=\"color:purple\">[Q]: What data/information do we need, in order to perform trio-phasing?</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[A]:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span style=\"color:purple\">[Q]: In order to run PBT, we need to supply it a \".ped\" file, as described here: http://www.gwaspi.org/?page_id=145 , where at least the first 6 fields are filled in. Why do we need this file for PBT? </span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[A]: "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span style=\"color:purple\">[Q]: How does PBT mark phased genotypes?</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[A]: "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A VCF file with all genotypes across brca1 and brca2 for the na12878 trio as well as an appropriate ped file can be found in the datafiles/COO5/ folder. \n",
    "\n",
    "\n",
    "Run GATK PBT accordingly, in order to obtain phased genotypes for NA12878 across these 2 genes (i.e.: see how and specify the region of interest -- brca1 and brca2 genes -- to GATK so it does not attempt to perform the analysis genome-wide)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "-java -jar $GATK -T PhaseByTransmission \\\n",
    "    -R  \\\n",
    "    -ped  \\\n",
    "    -V  \\\n",
    "    -L  \\\n",
    "    -o "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span style=\"color:purple\">[Q]: How many genotypes were NOT phased across the two genes, for NA12878, and why? (if there are any)</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[A]:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.Read-backed phasing of genetic variants (using GATK ReadBackedPhasing - RBP)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Look at the GATK RBP documentation for more details on what it needs and what it does, here: \n",
    "\n",
    "https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_gatk_tools_walkers_phasing_ReadBackedPhasing.php"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span style=\"color:purple\">[Q]: What data/information do we need in order to obtain read-back phased genotypes for an individual ?</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[A]: "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span style=\"color:purple\">[Q]: How does RBP mark a phased genotype ? </span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[A]:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run ReadBackedPhasing on the same set of variants that you ran PBT on, for sample NA12878, and look at output (use the whole genome aligned read-data that you aligned during the first exercise sessions):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "-java -jar $GATK -T ReadBackedPhasing \\\n",
    "    -R  \\\n",
    "    -sampleToPhase  \\\n",
    "    -V  \\\n",
    "    -I  \\\n",
    "    -L  \\\n",
    "    -o "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"grep\", as used below (ie.: with -v argument) retrieves all lines from a file (for our purpose, the vcf that you create by running RBP) that do NOT contain a certain TEXT. Replace TEXT in the command below, such that you retrieve all the RBP sites that are NOT phased. The subsequent \"wc -l\" command will count the number of such lines, giving you the final number of unphased genotypes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "grep -E -v '^#' | grep -v TEXT | wc -l"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span style=\"color:purple\">[Q]: What tool phased more variant sites in our case, PBT or RBP ? Which one would you expect to phase more sites and why? </span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[A]:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because RBP phases variant sites by overlapping reads that span adjacent variants, it will not offer global phasing (like PBT), but rather local phasing, resulting into local haplotypes, where one knows (or can infer from the output) how the alleles of 3 variants are combined on haplotypes and how the allleles of 10 other variants are combined onto haplotypes and so on. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span style=\"color:purple\">[Q]: (may require a piece of paper and a pen!) Extract the first three (independent) haplotype streches from the RBP output and represent these haplotypes as the way PBT does (i.e.: A|B). Then look at the PBT phasing for each of these 3 first haplotype stretches and evaluate whether the two phasings match or not  </span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[A]:"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Bash",
   "language": "bash",
   "name": "bash"
  },
  "language_info": {
   "codemirror_mode": "shell",
   "file_extension": ".sh",
   "mimetype": "text/x-sh",
   "name": "bash"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}