example-files/notebooks/aws-1000genomes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exploratory Data Analysis of Genomic Datasets with ADAM and Mango"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Configuring ADAM and Mango on EMR"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mango uses docker containers to be run easily on EMR. To get everything setup and installed, follow EMR documentation at https://bdg-mango.readthedocs.io/en/latest/cloud/emr.html."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Loading Data from the 1000 Genomes Project\n",
    "\n",
    "In this tutorial, we will use ADAM and Mango to discover interesting variants in the child of a 1000 Genomes trio.\n",
    "\n",
    "First, let’s import ADAM and Mango modules, as well as any Spark modules we need:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import ADAM modules\n",
    "from bdgenomics.adam.adamContext import ADAMContext\n",
    "from bdgenomics.adam.ds import AlignmentDataset, CoverageDataset\n",
    "from bdgenomics.adam.stringency import LENIENT, _toJava\n",
    "\n",
    "# Import Mango modules\n",
    "from bdgenomics.mango.alignments import *\n",
    "from bdgenomics.mango.coverage import CoverageDistribution\n",
    "\n",
    "# Import Spark modules\n",
    "from pyspark.sql import functions as sf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we will create an ADAMContext. ADAMContext allows us to load and manipulate genomic data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create ADAM Context\n",
    "ac = ADAMContext(spark)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Variant Analysis with Spark SQL"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this analysis, we will view a trio (NA19685, NA19661, and NA19660) and search for variants that are present in the child but not present in the parents. These are interesting regions, as they may indicate sights of de novo variation that may contribute to multiple disorders.\n",
    "\n",
    "First, we will load in a subset of variant data from chromosome 17:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pathPrefix = 's3://1000genomes/phase1/analysis_results/integrated_call_sets/'\n",
    "\n",
    "genotypesPath = pathPrefix + 'ALL.chr17.integrated_phase1_v3.20101123.snps_indels_svs.genotypes.vcf.gz'\n",
    "genotypes = ac.loadGenotypes(genotypesPath)\n",
    "\n",
    "genotypes_df  = genotypes.toDF()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can take a look at the schema by printing the columns in the dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cache genotypes and show the schema\n",
    "genotypes_df.columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This genotypes dataset contains all samples from the 1000 Genomes Project. Therefore, we will next filter genotypes to only consider samples that are in the NA19685 trio."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# trio IDs\n",
    "IDs = ['NA19685','NA19661','NA19660']\n",
    "\n",
    "trio_df = genotypes_df.filter(genotypes_df[\"sampleId\"].isin(IDs))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will next add a new column to our dataframe that determines the genomic location of each variant. This is defined by the chromosome (referenceName) and the start and end position of the variant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add ReferenceRegion column and group by referenceRegion\n",
    "trios_with_referenceRegion = trio_df.withColumn('ReferenceRegion', \n",
    "                    sf.concat(sf.col('referenceName'),sf.lit(':'), sf.col('start'), sf.lit('-'), sf.col('end')))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we want to query our dataset to find de novo variants. But first, we must register our dataframe with Spark SQL."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#  Register df with Spark SQL\n",
    "trios_with_referenceRegion.createOrReplaceTempView(\"trios\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that our dataframe is registered, we can run SQL queries on it. For our first query, we will select the names of a subset of variants belonging to sample NA19685 that have at least one alternative (ALT) allele."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# filter by alleles. This is a list of variant names that have an alternate allele for the child\n",
    "alternate_variant_sites = spark.sql(\"SELECT variant.names[0] AS snp FROM trios \\\n",
    "                                    WHERE array_contains(alleles, 'ALT') AND sampleId == 'NA19685'\") \n",
    "\n",
    "collected_sites = list(map(lambda x: x.snp, alternate_variant_sites.take(100)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For our next query, we will filter a subset of sites in which the parents have both reference alleles. We then filter these variants by the set produced above from the child."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get parent records and filter by only REF locations for variant names that were found in the child with an ALT\n",
    "filtered1 = spark.sql(\"SELECT * FROM trios WHERE sampleId == 'NA19661' or sampleId == 'NA19660' \\\n",
    "            AND !array_contains(alleles, 'ALT')\")\n",
    "\n",
    "\n",
    "filtered2 = filtered1.filter(filtered1[\"variant.names\"][0].isin(collected_sites))\n",
    "\n",
    "snp_counts = filtered2.take(100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# view snp names as a list\n",
    "set([x.variant.names[0] for x in snp_counts])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Working with Alignment Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we can explore these specific variant sites in the raw genomic alignment data. First, let’s load in the data for the NA19685 trio:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load in NA19685 exome from s3a\n",
    "childReadsPath = 's3a://1000genomes/phase1/data/NA19685/exome_alignment/NA19685.mapped.illumina.mosaik.MXL.exome.20110411.bam'\n",
    "parent1ReadsPath = 's3a://1000genomes/phase1/data/NA19660/exome_alignment/NA19660.mapped.illumina.mosaik.MXL.exome.20110411.bam'\n",
    "parent2ReadsPath = 's3a://1000genomes/phase1/data/NA19661/exome_alignment/NA19661.mapped.illumina.mosaik.MXL.exome.20110411.bam'\n",
    "\n",
    "childReads = ac.loadAlignments(childReadsPath, stringency=LENIENT)\n",
    "parent1Reads = ac.loadAlignments(parent1ReadsPath, stringency=LENIENT)\n",
    "parent2Reads = ac.loadAlignments(parent2ReadsPath, stringency=LENIENT)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Quality Control of Alignment Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One popular analysis to visually re-affirm the quality of genomic alignment data is by viewing coverage distribution. Coverage distribution gives us an idea of the read coverage we have across a sample. Next, we will generate a sample coverage distribution plot for the child alignment data on chromosome 17."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate read coverage\n",
    "# Takes 2-3 minutes\n",
    "childCoverage = childReads.transform(lambda x: x.filter(x.referenceName == \"17\")).toCoverage()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that coverage data is calculated and cached, we will compute the coverage distribution of all three samples and plot the coverage distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate coverage distribution\n",
    "# You can check the progress in the SparkUI by navigating to \n",
    "# <PUBLIC_MASTER_DNS>:8088 and clicking on the currently running Spark application.\n",
    "cd = CoverageDistribution(spark, childCoverage, bin_size = 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax, results = cd.plotDistributions(normalize=True, cumulative=False)\n",
    "\n",
    "ax.set_title(\"Coverage Distribution\")\n",
    "ax.set_ylabel(\"Counts\")\n",
    "ax.set_xlabel(\"Coverage Depth\")\n",
    "ax.set_xscale(\"log\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we are done with coverage, we can unpersist these datasets to clear space in memory for the next analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "childCoverage.unpersist()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Viewing Sites with Missense Variants in the Proband\n",
    "\n",
    "After verifying alignment data and filtering variants, we have 4 genes with potential missense mutations in the proband, including YBX2, ZNF286B, KSR1, and GNA13. We can visually verify these sites by filtering and viewing the raw reads of the child and parents.\n",
    "\n",
    "First, let's view the child reads. If we zoom in to the location of the GNA13 variant (63052580-63052581) we can see a heterozygous T to A call."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# view missense variant at GNA13: 63052580-63052581 (SNP rs201316886) in child\n",
    "# Takes about 2 minutes to collect data from workers\n",
    "childViz = AlignmentSummary(spark, ac, childReads)\n",
    "contig = \"17\"\n",
    "start = 63052180\n",
    "end = 63052981\n",
    "\n",
    "childViz.viewPileup(contig, start, end)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like there indeed is a variant at this position, possibly a heterozygous SNP with alternate allele A. Let's look at the parent data to verify this variant does not appear in the parents."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# view missense variant at GNA13: 63052580-63052581 in parent 1\n",
    "parent1Viz = AlignmentSummary(spark, ac, parent1Reads)\n",
    "contig = \"17\"\n",
    "start = 63052180\n",
    "end = 63052981\n",
    "\n",
    "parent1Viz.viewPileup(contig, start, end)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# view missense variant at GNA13: 63052580-63052581 in parent 2 \n",
    "parent2Viz = AlignmentSummary(spark, ac, parent2Reads)\n",
    "contig = \"17\"\n",
    "start = 63052180\n",
    "end = 63052981\n",
    "\n",
    "parent2Viz.viewPileup(contig, start, end)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This confirms our filter that this variant is indeed only present in the proband, but not the parents."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "To summarize, this post demonstrated how to setup and run ADAM and Mango in EMR. We demonstrated how to use these tools in an interactive notebook environment to explore the 1000 Genomes Dataset, a publicly available dataset on Amazon S3. We used these tools inspect 1000 Genomes data quality, query for interesting variants in the genome and validate results through visualization of raw datsets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}