03_glow_vcf_explore_merged.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import pyspark\n",
    "import glow"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "builder = pyspark.sql.SparkSession.builder.appName(\"GlowVCFExplore\") \\\n",
    "    .config(\"spark.jars.packages\", \"io.projectglow:glow-spark3_2.12:1.2.1\") \\\n",
    "    .config(\"spark.hadoop.io.compression.codecs\", \"io.projectglow.sql.util.BGZFCodec\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "spark = builder.getOrCreate()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "spark = glow.register(spark)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spark"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# Page Break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Explore bcBio both Somatic and Germline VCFs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "22,23d21\n",
      "<  |-- INFO_FREQ: array (nullable = true)\n",
      "<  |    |-- element: string (containsNull = true)\n",
      "37,38d34\n",
      "<  |-- INFO_SOMTYPE: array (nullable = true)\n",
      "<  |    |-- element: string (containsNull = true)\n"
     ]
    }
   ],
   "source": [
    "# Observe that somatic schema has 2 additional columns `INFO_FREQ` and `INFO_SOMTYPE`\n",
    "!diff schema_bcbio_giab_somatic.txt schema_bcbio_giab_germline.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# This time we just point to data directory\n",
    "bcbio_src = \"./data/bcbio_giab_somatic/*.vcf.gz\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bcbio_df = spark.read.format(\"vcf\").load(bcbio_src)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "root\n",
      " |-- contigName: string (nullable = true)\n",
      " |-- start: long (nullable = true)\n",
      " |-- end: long (nullable = true)\n",
      " |-- names: array (nullable = true)\n",
      " |    |-- element: string (containsNull = true)\n",
      " |-- referenceAllele: string (nullable = true)\n",
      " |-- alternateAlleles: array (nullable = true)\n",
      " |    |-- element: string (containsNull = true)\n",
      " |-- qual: double (nullable = true)\n",
      " |-- filters: array (nullable = true)\n",
      " |    |-- element: string (containsNull = true)\n",
      " |-- splitFromMultiAllelic: boolean (nullable = true)\n",
      " |-- INFO_platformnames: array (nullable = true)\n",
      " |    |-- element: string (containsNull = true)\n",
      " |-- INFO_callsetwithotheruniqgenopassing: array (nullable = true)\n",
      " |    |-- element: string (containsNull = true)\n",
      " |-- INFO_callsetnames: array (nullable = true)\n",
      " |    |-- element: string (containsNull = true)\n",
      " |-- INFO_AC: array (nullable = true)\n",
      " |    |-- element: integer (containsNull = true)\n",
      " |-- INFO_FREQ: array (nullable = true)\n",
      " |    |-- element: string (containsNull = true)\n",
      " |-- INFO_varType: string (nullable = true)\n",
      " |-- INFO_DPSum: integer (nullable = true)\n",
      " |-- INFO_datasetsmissingcall: array (nullable = true)\n",
      " |    |-- element: string (containsNull = true)\n",
      " |-- INFO_AN: integer (nullable = true)\n",
      " |-- INFO_callsets: integer (nullable = true)\n",
      " |-- INFO_callsetwiththisuniqgenopassing: array (nullable = true)\n",
      " |    |-- element: string (containsNull = true)\n",
      " |-- INFO_callable: array (nullable = true)\n",
      " |    |-- element: string (containsNull = true)\n",
      " |-- INFO_difficultregion: array (nullable = true)\n",
      " |    |-- element: string (containsNull = true)\n",
      " |-- INFO_datasets: integer (nullable = true)\n",
      " |-- INFO_SOMTYPE: array (nullable = true)\n",
      " |    |-- element: string (containsNull = true)\n",
      " |-- INFO_platformbias: array (nullable = true)\n",
      " |    |-- element: string (containsNull = true)\n",
      " |-- INFO_platforms: integer (nullable = true)\n",
      " |-- INFO_arbitrated: string (nullable = true)\n",
      " |-- INFO_filt: array (nullable = true)\n",
      " |    |-- element: string (containsNull = true)\n",
      " |-- INFO_datasetnames: array (nullable = true)\n",
      " |    |-- element: string (containsNull = true)\n",
      " |-- genotypes: array (nullable = true)\n",
      " |    |-- element: struct (containsNull = true)\n",
      " |    |    |-- sampleId: string (nullable = true)\n",
      " |    |    |-- conditionalQuality: integer (nullable = true)\n",
      " |    |    |-- alleleDepths: array (nullable = true)\n",
      " |    |    |    |-- element: integer (containsNull = true)\n",
      " |    |    |-- phased: boolean (nullable = true)\n",
      " |    |    |-- calls: array (nullable = true)\n",
      " |    |    |    |-- element: integer (containsNull = true)\n",
      " |    |    |-- depth: integer (nullable = true)\n",
      " |    |    |-- IGT: string (nullable = true)\n",
      " |    |    |-- IPS: string (nullable = true)\n",
      " |    |    |-- ADALL: array (nullable = true)\n",
      " |    |    |    |-- element: integer (containsNull = true)\n",
      " |    |    |-- PS: integer (nullable = true)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Observe that Spark has merged the schema i.e. Union of columns from both VCF headers\n",
    "bcbio_df.printSchema()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bcbio_df.createOrReplaceTempView(\"vcf_table\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+--------------------+--------------------+-------+\n",
      "|            col_name|           data_type|comment|\n",
      "+--------------------+--------------------+-------+\n",
      "|          contigName|              string|   null|\n",
      "|               start|              bigint|   null|\n",
      "|                 end|              bigint|   null|\n",
      "|               names|       array<string>|   null|\n",
      "|     referenceAllele|              string|   null|\n",
      "|    alternateAlleles|       array<string>|   null|\n",
      "|                qual|              double|   null|\n",
      "|             filters|       array<string>|   null|\n",
      "|splitFromMultiAll...|             boolean|   null|\n",
      "|  INFO_platformnames|       array<string>|   null|\n",
      "|INFO_callsetwitho...|       array<string>|   null|\n",
      "|   INFO_callsetnames|       array<string>|   null|\n",
      "|             INFO_AC|          array<int>|   null|\n",
      "|           INFO_FREQ|       array<string>|   null|\n",
      "|        INFO_varType|              string|   null|\n",
      "|          INFO_DPSum|                 int|   null|\n",
      "|INFO_datasetsmiss...|       array<string>|   null|\n",
      "|             INFO_AN|                 int|   null|\n",
      "|       INFO_callsets|                 int|   null|\n",
      "|INFO_callsetwitht...|       array<string>|   null|\n",
      "|       INFO_callable|       array<string>|   null|\n",
      "|INFO_difficultregion|       array<string>|   null|\n",
      "|       INFO_datasets|                 int|   null|\n",
      "|        INFO_SOMTYPE|       array<string>|   null|\n",
      "|   INFO_platformbias|       array<string>|   null|\n",
      "|      INFO_platforms|                 int|   null|\n",
      "|     INFO_arbitrated|              string|   null|\n",
      "|           INFO_filt|       array<string>|   null|\n",
      "|   INFO_datasetnames|       array<string>|   null|\n",
      "|           genotypes|array<struct<samp...|   null|\n",
      "+--------------------+--------------------+-------+\n",
      "\n"
     ]
    }
   ],
   "source": [
    "spark.sql(\"describe vcf_table\").show(n=1000, truncate=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING\t2022-10-07 23:28:40\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "WARNING\t2022-10-07 23:28:40\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "WARNING\t2022-10-07 23:28:40\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "WARNING\t2022-10-07 23:28:40\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "WARNING\t2022-10-07 23:28:43\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "[Stage 2:============================================>              (3 + 1) / 4]\r"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+------------------+\n",
      "|number_of_variants|\n",
      "+------------------+\n",
      "|           3206719|\n",
      "+------------------+\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "                                                                                \r"
     ]
    }
   ],
   "source": [
    "# Observe that total variants count from both VCFs in single table\n",
    "spark.sql(\"select count(1) as number_of_variants from vcf_table\").show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING\t2022-10-07 23:28:49\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "WARNING\t2022-10-07 23:28:49\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "WARNING\t2022-10-07 23:28:49\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "WARNING\t2022-10-07 23:28:49\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "WARNING\t2022-10-07 23:28:52\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "[Stage 5:============================================>              (3 + 1) / 4]\r"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+----------+\n",
      "|contigName|\n",
      "+----------+\n",
      "|      chr1|\n",
      "|     chr10|\n",
      "|     chr11|\n",
      "|     chr12|\n",
      "|     chr13|\n",
      "|     chr14|\n",
      "|     chr15|\n",
      "|     chr16|\n",
      "|     chr17|\n",
      "|     chr18|\n",
      "|     chr19|\n",
      "|      chr2|\n",
      "|     chr20|\n",
      "|     chr21|\n",
      "|     chr22|\n",
      "|      chr3|\n",
      "|      chr4|\n",
      "|      chr5|\n",
      "|      chr6|\n",
      "|      chr7|\n",
      "|      chr8|\n",
      "|      chr9|\n",
      "+----------+\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "                                                                                \r"
     ]
    }
   ],
   "source": [
    "spark.sql(\"select distinct contigName from vcf_table order by contigName\").show(46)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+----------+------+------+\n",
      "|contigName| start|   end|\n",
      "+----------+------+------+\n",
      "|      chr1|817185|817186|\n",
      "|      chr1|817340|817341|\n",
      "|      chr1|817888|817889|\n",
      "|      chr1|818801|818802|\n",
      "|      chr1|818811|818812|\n",
      "|      chr1|818953|818954|\n",
      "|      chr1|819122|819123|\n",
      "|      chr1|819583|819584|\n",
      "|      chr1|824319|824320|\n",
      "|      chr1|824456|824457|\n",
      "|      chr1|825531|825532|\n",
      "|      chr1|825766|825767|\n",
      "|      chr1|826576|826577|\n",
      "|      chr1|826892|826893|\n",
      "|      chr1|827208|827209|\n",
      "|      chr1|827211|827212|\n",
      "|      chr1|827220|827221|\n",
      "|      chr1|827251|827252|\n",
      "|      chr1|828013|828014|\n",
      "|      chr1|830724|830725|\n",
      "+----------+------+------+\n",
      "only showing top 20 rows\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING\t2022-10-07 23:28:58\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n"
     ]
    }
   ],
   "source": [
    "spark.sql(\"select contigName, start, end from vcf_table\").show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING\t2022-10-07 23:29:02\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "WARNING\t2022-10-07 23:29:02\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "WARNING\t2022-10-07 23:29:02\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "WARNING\t2022-10-07 23:29:02\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "WARNING\t2022-10-07 23:29:04\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "[Stage 9:============================================>              (3 + 1) / 4]\r"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+----------+----------+\n",
      "|contigName|num_of_pos|\n",
      "+----------+----------+\n",
      "|      chr2|    271356|\n",
      "|      chr1|    267876|\n",
      "|      chr3|    243088|\n",
      "|      chr6|    220641|\n",
      "|      chr4|    213774|\n",
      "|      chr5|    204282|\n",
      "|      chr7|    186041|\n",
      "|     chr11|    172760|\n",
      "|     chr10|    169586|\n",
      "|      chr8|    165373|\n",
      "|     chr12|    154561|\n",
      "|      chr9|    147990|\n",
      "|     chr13|    137088|\n",
      "|     chr14|    112042|\n",
      "|     chr15|     96779|\n",
      "|     chr17|     83405|\n",
      "|     chr18|     80893|\n",
      "|     chr20|     71827|\n",
      "|     chr19|     61776|\n",
      "|     chr16|     60843|\n",
      "|     chr21|     47651|\n",
      "|     chr22|     37087|\n",
      "+----------+----------+\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "                                                                                \r"
     ]
    }
   ],
   "source": [
    "spark.sql(\"select contigName, count(end) as num_of_pos from vcf_table group by contigName order by num_of_pos desc\").show(46)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+---------------+----------------+----------------------------+\n",
      "|referenceAllele|alternateAlleles|array_size(alternateAlleles)|\n",
      "+---------------+----------------+----------------------------+\n",
      "|              G|             [A]|                           1|\n",
      "|              A|             [G]|                           1|\n",
      "|              C|             [G]|                           1|\n",
      "|              A|             [G]|                           1|\n",
      "|              A|             [G]|                           1|\n",
      "|              T|             [C]|                           1|\n",
      "|              G|             [A]|                           1|\n",
      "|              C|             [T]|                           1|\n",
      "|              T|             [C]|                           1|\n",
      "|              T|             [A]|                           1|\n",
      "|              C|             [T]|                           1|\n",
      "|              T|             [C]|                           1|\n",
      "|              A|            [AT]|                           1|\n",
      "|              G|             [A]|                           1|\n",
      "|              G|             [C]|                           1|\n",
      "|              C|             [G]|                           1|\n",
      "|              T|             [C]|                           1|\n",
      "|              T|             [A]|                           1|\n",
      "|              G|             [A]|                           1|\n",
      "|              T|             [A]|                           1|\n",
      "+---------------+----------------+----------------------------+\n",
      "only showing top 20 rows\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING\t2022-10-07 23:29:10\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n"
     ]
    }
   ],
   "source": [
    "spark.sql(\"select referenceAllele, alternateAlleles, array_size(alternateAlleles) from vcf_table\").show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING\t2022-10-07 23:29:14\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "WARNING\t2022-10-07 23:29:14\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "WARNING\t2022-10-07 23:29:14\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "WARNING\t2022-10-07 23:29:14\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "WARNING\t2022-10-07 23:29:16\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "[Stage 13:==============>                                           (1 + 3) / 4]\r"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+---------------+----------------+-----------+\n",
      "|referenceAllele|alternateAlleles|num_of_snps|\n",
      "+---------------+----------------+-----------+\n",
      "|              C|             [T]|     484436|\n",
      "|              G|             [A]|     484020|\n",
      "|              T|             [C]|     463525|\n",
      "|              A|             [G]|     462312|\n",
      "|              G|             [C]|     119203|\n",
      "|              C|             [G]|     118972|\n",
      "|              G|             [T]|     117604|\n",
      "|              C|             [A]|     117027|\n",
      "|              A|             [C]|     114401|\n",
      "|              T|             [G]|     113450|\n",
      "|              T|             [A]|      96640|\n",
      "|              A|             [T]|      96202|\n",
      "+---------------+----------------+-----------+\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "                                                                                \r"
     ]
    }
   ],
   "source": [
    "spark.sql(\"select referenceAllele, alternateAlleles, count(*) as num_of_snps \\\n",
    "from vcf_table \\\n",
    "where \\\n",
    "    char_length(referenceAllele) = 1 and \\\n",
    "    array_size(alternateAlleles) = 1 and \\\n",
    "    char_length(alternateAlleles[0]) = 1 \\\n",
    "    group by referenceAllele, alternateAlleles \\\n",
    "    order by num_of_snps desc\").show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+----------+------+------+---------------+----------------+---------+------------+\n",
      "|contigName|start |end   |referenceAllele|alternateAlleles|sampleId |alleleDepths|\n",
      "+----------+------+------+---------------+----------------+---------+------------+\n",
      "|chr1      |817185|817186|G              |[A]             |[NA12878]|[[78, 454]] |\n",
      "|chr1      |817340|817341|A              |[G]             |[NA12878]|[[107, 342]]|\n",
      "|chr1      |817888|817889|C              |[G]             |[NA12878]|[[74, 220]] |\n",
      "|chr1      |818801|818802|A              |[G]             |[NA12878]|[[0, 202]]  |\n",
      "|chr1      |818811|818812|A              |[G]             |[NA12878]|[[0, 190]]  |\n",
      "|chr1      |818953|818954|T              |[C]             |[NA12878]|[[0, 246]]  |\n",
      "|chr1      |819122|819123|G              |[A]             |[NA12878]|[[110, 363]]|\n",
      "|chr1      |819583|819584|C              |[T]             |[NA12878]|[[91, 378]] |\n",
      "|chr1      |824319|824320|T              |[C]             |[NA12878]|[[92, 391]] |\n",
      "|chr1      |824456|824457|T              |[A]             |[NA12878]|[[0, 355]]  |\n",
      "|chr1      |825531|825532|C              |[T]             |[NA12878]|[[0, 243]]  |\n",
      "|chr1      |825766|825767|T              |[C]             |[NA12878]|[[0, 244]]  |\n",
      "|chr1      |826576|826577|A              |[AT]            |[NA12878]|[[0, 0]]    |\n",
      "|chr1      |826892|826893|G              |[A]             |[NA12878]|[[60, 224]] |\n",
      "|chr1      |827208|827209|G              |[C]             |[NA12878]|[[0, 259]]  |\n",
      "|chr1      |827211|827212|C              |[G]             |[NA12878]|[[0, 263]]  |\n",
      "|chr1      |827220|827221|T              |[C]             |[NA12878]|[[0, 259]]  |\n",
      "|chr1      |827251|827252|T              |[A]             |[NA12878]|[[0, 257]]  |\n",
      "|chr1      |828013|828014|G              |[A]             |[NA12878]|[[0, 0]]    |\n",
      "|chr1      |830724|830725|T              |[A]             |[NA12878]|[[0, 171]]  |\n",
      "+----------+------+------+---------------+----------------+---------+------------+\n",
      "only showing top 20 rows\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING\t2022-10-07 23:29:26\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n"
     ]
    }
   ],
   "source": [
    "spark.sql(\"select contigName, start, end, referenceAllele, alternateAlleles, genotypes.sampleId, genotypes.alleleDepths \\\n",
    "from vcf_table\").show(truncate=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING\t2022-10-07 23:29:29\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "WARNING\t2022-10-07 23:29:32\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "WARNING\t2022-10-07 23:29:32\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "WARNING\t2022-10-07 23:29:32\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "WARNING\t2022-10-07 23:29:33\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "[Stage 18:===================>                                      (1 + 2) / 3]\r"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+----------+------+------+---------------+----------------+---------+\n",
      "|contigName|start |end   |referenceAllele|alternateAlleles|sampleId |\n",
      "+----------+------+------+---------------+----------------+---------+\n",
      "|chr1      |817185|817186|G              |[A]             |[NA12878]|\n",
      "+----------+------+------+---------------+----------------+---------+\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "                                                                                \r"
     ]
    }
   ],
   "source": [
    "spark.sql(\"select contigName, start, end, referenceAllele, alternateAlleles, genotypes.sampleId from vcf_table \\\n",
    "where contigName = 'chr1' and end = 817186\").show(truncate=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+----------+------+------+---------------+----------------+---------------------+\n",
      "|contigName|start |end   |referenceAllele|alternateAlleles|genotypes.sampleId[0]|\n",
      "+----------+------+------+---------------+----------------+---------------------+\n",
      "|chr1      |817185|817186|G              |[A]             |NA12878              |\n",
      "|chr1      |817340|817341|A              |[G]             |NA12878              |\n",
      "|chr1      |817888|817889|C              |[G]             |NA12878              |\n",
      "|chr1      |818801|818802|A              |[G]             |NA12878              |\n",
      "|chr1      |818811|818812|A              |[G]             |NA12878              |\n",
      "|chr1      |818953|818954|T              |[C]             |NA12878              |\n",
      "|chr1      |819122|819123|G              |[A]             |NA12878              |\n",
      "|chr1      |819583|819584|C              |[T]             |NA12878              |\n",
      "|chr1      |824319|824320|T              |[C]             |NA12878              |\n",
      "|chr1      |824456|824457|T              |[A]             |NA12878              |\n",
      "|chr1      |825531|825532|C              |[T]             |NA12878              |\n",
      "|chr1      |825766|825767|T              |[C]             |NA12878              |\n",
      "|chr1      |826576|826577|A              |[AT]            |NA12878              |\n",
      "|chr1      |826892|826893|G              |[A]             |NA12878              |\n",
      "|chr1      |827208|827209|G              |[C]             |NA12878              |\n",
      "|chr1      |827211|827212|C              |[G]             |NA12878              |\n",
      "|chr1      |827220|827221|T              |[C]             |NA12878              |\n",
      "|chr1      |827251|827252|T              |[A]             |NA12878              |\n",
      "|chr1      |828013|828014|G              |[A]             |NA12878              |\n",
      "|chr1      |830724|830725|T              |[A]             |NA12878              |\n",
      "+----------+------+------+---------------+----------------+---------------------+\n",
      "only showing top 20 rows\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING\t2022-10-07 23:29:38\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n"
     ]
    }
   ],
   "source": [
    "spark.sql(\"select contigName, start, end, referenceAllele, alternateAlleles, genotypes.sampleId[0] from vcf_table \\\n",
    "where genotypes.sampleId[0] = 'NA12878'\").show(truncate=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# Page Break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Filtering Somatic or Germline\n",
    "\n",
    "* Single merged `vcf_table` is great!\n",
    "* However. Can I still filter, say, I wanted to query Somatic records only or vice versa?\n",
    "* Recall that `INFO_SOMTYPE` or `INFO_FREQ` columns only present in Somatic VCF header.\n",
    "* Hence, we can approximate `NULL` record present of either column as pivotal data filter.\n",
    "* Let try example with `INFO_SOMTYPE` column."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING\t2022-10-07 23:29:46\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "WARNING\t2022-10-07 23:29:46\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "WARNING\t2022-10-07 23:29:46\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "WARNING\t2022-10-07 23:29:46\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "WARNING\t2022-10-07 23:29:48\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "[Stage 20:==============>                                           (1 + 3) / 4]\r"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+------------------+\n",
      "|number_of_variants|\n",
      "+------------------+\n",
      "|           3206719|\n",
      "+------------------+\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "                                                                                \r"
     ]
    }
   ],
   "source": [
    "# Count total records\n",
    "spark.sql(\"select count(1) as number_of_variants from vcf_table\").show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING\t2022-10-07 23:29:53\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "WARNING\t2022-10-07 23:29:53\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "WARNING\t2022-10-07 23:29:53\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "WARNING\t2022-10-07 23:29:53\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "WARNING\t2022-10-07 23:29:55\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "[Stage 23:=============================>                            (2 + 2) / 4]\r"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+--------------------------+\n",
      "|number_of_variants_somatic|\n",
      "+--------------------------+\n",
      "|                   1082945|\n",
      "+--------------------------+\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "                                                                                \r"
     ]
    }
   ],
   "source": [
    "# Filter Somatic only records\n",
    "spark.sql(\"select count(1) as number_of_variants_somatic from vcf_table where INFO_SOMTYPE is not null\").show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING\t2022-10-07 23:30:01\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "WARNING\t2022-10-07 23:30:01\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "WARNING\t2022-10-07 23:30:01\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "WARNING\t2022-10-07 23:30:01\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "WARNING\t2022-10-07 23:30:03\tAsciiLineReader\tCreating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream\n",
      "[Stage 26:=============================>                            (2 + 2) / 4]\r"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+---------------------------+\n",
      "|number_of_variants_germline|\n",
      "+---------------------------+\n",
      "|                    2123774|\n",
      "+---------------------------+\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "                                                                                \r"
     ]
    }
   ],
   "source": [
    "# Filter Germline only records\n",
    "spark.sql(\"select count(1) as number_of_variants_germline from vcf_table where INFO_SOMTYPE is null\").show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Summary\n",
    "\n",
    "* We can process VCF by their study type Somatic or Germline calling or both.\n",
    "* Assumption is that within each \"best practice\" BioInformatics Pipeline; it should generate similar VCF structure with minor differences in header annotation.\n",
    "* Couple of strategy possible:\n",
    "    * by arranging all Somatic VCF type of the same BioInfo Pipeline output into one table\n",
    "    * similarly, all Germline VCF type of the same BioInfo Pipeline output into one table\n",
    "    * if we merged VCF, make sure to have very discriminator column that can filter data records better\n",
    "        * depends on data pipline setup, this discriminator column can be inserted during post-processing VCF files\n",
    "        * or, could add as part of BioInfo Pipeline VCF annotation process\n",
    "* We can prescribe Spark/Glow to auto-discover \"schema\" out of VCF. Hence, \"schema evolution\" is possible.\n",
    "* Or, we can prescribe \"pre-defined schema\" to Spark/Glow during parsing. Hence, enforcing \"strict schema\".\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# Page Break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Stop Spark Session"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "spark.stop()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# Continue to next notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}