sample-notebooks/genomics-platinum-genomes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Getting the Illumina Platinum Genomes from Azure Open Datasets and Doing Initial Analysis "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Jupyter notebooks are a great tool for data scientists who is working on Genomics data analysis. We will demonstrate Azure Jupyter notebook usage via GATK and Picard with Azure Open Datasets. \n",
    "\n",
    "**Here is the coverage of this notebook:**\n",
    "\n",
    "1. Annotate genotypes using VariantFiltration\n",
    "2. Select Specific Variants\n",
    "3. Filter the relevant variants- no calls OR specific regions\n",
    "4. Perform concordance analysis\n",
    "5. Convert the final VCF files to a table \n",
    "\n",
    "**Dependencies:**\n",
    "\n",
    "This notebook requires the following libraries:\n",
    "\n",
    "- Azure storage `pip install azure-storage-blob`\n",
    "\n",
    "- numpy `pip install numpy`\n",
    "\n",
    "- Genome Analysis Toolkit (GATK) (*Users need to download GATK from Broad Institute's webpage into the same compute environment with this notebook: https://github.com/broadinstitute/gatk/releases*)\n",
    "\n",
    "**Important information: This notebook is using Python 3.6 kernel**\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Getting the Genomics data from Azure Open Datasets\n",
    "\n",
    "Several public genomics data has been uploaded as an Azure Open Dataset [here](https://azure.microsoft.com/services/open-datasets/catalog/). We create a blob service linked to this open datasets. You can find examples of data calling procedure from Azure Open Dataset for `Illumina Platinum Genomes` datasets in below:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Downloading the specific 'Illumina Platinum Genomes'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
    "source": [
      "import os\n",
      "import uuid\n",
      "import sys\n",
      "from azure.storage.blob import BlockBlobService, PublicAccess\n",
      "\n",
      "blob_service_client = BlockBlobService(account_name='datasetplatinumgenomes', sas_token='sv=2019-02-02&se=2050-01-01T08%3A00%3A00Z&si=prod&sr=c&sig=FFfZ0QaDcnEPQmWsshtpoYOjbzd4jtwIWeK%2Fc4i9MqM%3D')     \n",
      "blob_service_client.get_blob_to_path('dataset/2017-1.0/hg38/small_variants/NA12877', 'NA12877.vcf.gz', './NA12877.vcf.gz')"
    ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Annotate genotypes using VariantFiltration\n",
    "\n",
    "**Important note: Please check your GATK is running on your system.**\n",
    "\n",
    "If we want to filter heterozygous genotypes, we use VariantFiltration's `--genotype-filter-expression isHet == 1` option. We can specify the annotation value for the tool to label the heterozygous genotypes with with the `--genotype-filter-name` option. Here, this parameter's value is set to `isHetFilter`. In our first example, we used `NA12877.vcf.gz` from Illimina Platinum Genomes but users can use any vcf files from other datasets:`Platinum Genomes`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "run gatk VariantFiltration -V NA12877.vcf.gz -O outputannot.vcf --genotype-filter-expression \"isHet == 1\" --genotype-filter-name \"isHetFilter\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Select Specific Variants"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Select a subset of variants from a VCF file\n",
    "This tool makes it possible to select a subset of variants based on various criteria in order to facilitate certain analyses. Examples of such analyses include comparing and contrasting cases vs. controls, extracting variant or non-variant loci that meet certain requirements, or troubleshooting some unexpected results, to name a few.\n",
    "\n",
    "There are many different options for selecting subsets of variants from a larger call set:\n",
    "\n",
    "Extract one or more samples from a callset based on either a complete sample name or a pattern match.\n",
    "Specify criteria for inclusion that place thresholds on annotation values, **e.g. \"DP > 1000\" (depth of coverage greater than 1000x), \"AF < 0.25\" (sites with allele frequency less than 0.25)**.These criteria are written as \"JEXL expressions\", which are documented in the article about using JEXL expressions.\n",
    "Provide concordance or discordance tracks in order to include or exclude variants that are also present in other given callsets.\n",
    "Select variants based on criteria like their type (e.g. INDELs only), evidence of mendelian violation, filtering status, allelicity, etc.\n",
    "There are also several options for recording the original values of certain annotations which are recalculated when one subsets the new callset, trims alleles, etc.\n",
    "\n",
    "**Input**\n",
    "\n",
    "A variant call set in VCF format from which a subset can be selected.\n",
    "\n",
    "**Output**\n",
    "\n",
    "A new VCF file containing the selected subset of variants."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "run gatk SelectVariants -R Homo_sapiens_assembly38.fasta -V outputannot.vcf --select-type-to-include SNP --select-type-to-include INDEL -O selective.vcf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Transform filtered genotypes to no call \n",
    "\n",
    "Running SelectVariants with --set-filtered-gt-to-nocall will further transform the flagged genotypes with a null genotype call. \n",
    "\n",
    "This conversion is necessary because downstream tools do not parse the FORMAT-level filter field."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How can we filter the variants with with **'No call'**\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "run gatk SelectVariants -V outputannot.vcf --set-filtered-gt-to-nocall -O outputnocall.vcf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Check the Concordance of VCF file with Ground Truth"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Evaluate site-level concordance of an input VCF against a truth VCF.\n",
    "This tool evaluates two variant callsets against each other and produces a six-column summary metrics table. \n",
    "\n",
    "**This function will :**\n",
    "\n",
    "1. Stratifies SNP and INDEL calls\n",
    "2. Report true-positive,False-positive and false-negative calls\n",
    "3. Calculates sensitivity and precision\n",
    "\n",
    "The tool assumes all records in the --truth VCF are passing truth variants. For the -eval VCF, the tool uses only unfiltered passing calls.\n",
    "\n",
    "Optionally, the tool can be set to produce VCFs of the following variant records, annotated with each variant's concordance status:\n",
    "\n",
    "True positives and false negatives (i.e. all variants in the truth VCF): useful for calculating sensitivity\n",
    "\n",
    "True positives and false positives (i.e. all variants in the eval VCF): useful for obtaining a training data set for machine learning classifiers of artifacts\n",
    "\n",
    "**These output VCFs can be passed to VariantsToTable to produce a TSV file for statistical analysis in R or Python.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    " run gatk Concordance -R Homo_sapiens_assembly38.fasta -eval outputannot.vcf --truth outputnocall.vcf  --summary summary.tsv "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5. VariantsToTable"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extract fields from a VCF file to a tab-delimited table\n",
    "This tool extracts specified fields for each variant in a VCF file to a tab-delimited table, which may be easier to work with than a VCF. By default, the tool only extracts PASS or . (unfiltered) variants in the VCF file. Filtered variants may be included in the output by adding the --show-filtered flag. The tool can extract both INFO (i.e. site-level) fields and FORMAT (i.e. sample-level) fields.\n",
    "\n",
    "\n",
    "**INFO/site-level fields**\n",
    "\n",
    "Use the `-F` argument to extract INFO fields; each field will occupy a single column in the output file. The field can be any standard VCF column (e.g. CHROM, ID, QUAL) or any annotation name in the INFO field (e.g. AC, AF). The tool also supports the following additional fields:\n",
    "\n",
    "EVENTLENGTH (length of the event)\n",
    "TRANSITION (1 for a bi-allelic transition (SNP), 0 for bi-allelic transversion (SNP), -1 for INDELs and multi-allelics)\n",
    "HET (count of het genotypes)\n",
    "HOM-REF (count of homozygous reference genotypes)\n",
    "HOM-VAR (count of homozygous variant genotypes)\n",
    "NO-CALL (count of no-call genotypes)\n",
    "TYPE (type of variant, possible values are NO_VARIATION, SNP, MNP, INDEL, SYMBOLIC, and MIXED\n",
    "VAR (count of non-reference genotypes)\n",
    "NSAMPLES (number of samples)\n",
    "NCALLED (number of called samples)\n",
    "MULTI-ALLELIC (is this variant multi-allelic? true/false)\n",
    "\n",
    "\n",
    "**FORMAT/sample-level fields**\n",
    "\n",
    "Use the `-GF` argument to extract FORMAT/sample-level fields. The tool will create a new column per sample with the name \"SAMPLE_NAME.FORMAT_FIELD_NAME\" e.g. NA12877.GQ, NA12878.GQ.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Input**\n",
    "\n",
    "A VCF file to convert to a table\n",
    "\n",
    "**Output**\n",
    "\n",
    "A tab-delimited file containing the values of the requested fields in the VCF file.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "run gatk VariantsToTable -V NA12877.vcf.gz -F CHROM -F POS -F TYPE -F AC -F AD -F AF -GF DP -GF AD -O outputtable.table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# References\n",
    "\n",
    "1. VariantFiltration: https://gatk.broadinstitute.org/hc/en-us/articles/360036827111-VariantFiltration \n",
    "2. Select Variants:https://gatk.broadinstitute.org/hc/en-us/articles/360037052272-SelectVariants\n",
    "3. Concordance: https://gatk.broadinstitute.org/hc/en-us/articles/360041851651-Concordance\n",
    "4. Variants to table: https://gatk.broadinstitute.org/hc/en-us/articles/360036882811-VariantsToTable \n",
    "5. Illumina Platinum Genomes:https://www.illumina.com/platinumgenomes.html \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**END OF NOTEBOOK**"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "formats": "ipynb,md"
  },
  "kernelspec": {
   "display_name": "Python 3.6 - AzureML",
   "language": "python",
   "name": "python3-azureml"
  },
  "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}