03_glow_vcf_explore_merged_databricks.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "application/vnd.databricks.v1+cell": {
     "cellMetadata": {},
     "inputWidgets": {},
     "nuid": "7f3aeaa1-bef0-4c50-92c3-273d31ec84fa",
     "showTitle": false,
     "title": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# import pyspark\n",
    "# import glow"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "application/vnd.databricks.v1+cell": {
     "cellMetadata": {},
     "inputWidgets": {},
     "nuid": "1cade5e5-0621-4936-9d97-2dea69819e3b",
     "showTitle": false,
     "title": ""
    },
    "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": 0,
   "metadata": {
    "application/vnd.databricks.v1+cell": {
     "cellMetadata": {},
     "inputWidgets": {},
     "nuid": "97f07aeb-c647-4d02-b662-e38eee0ced37",
     "showTitle": false,
     "title": ""
    }
   },
   "outputs": [],
   "source": [
    "# spark = builder.getOrCreate()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "application/vnd.databricks.v1+cell": {
     "cellMetadata": {},
     "inputWidgets": {},
     "nuid": "39e7c162-9e0c-485c-85ac-3d819035cf97",
     "showTitle": false,
     "title": ""
    }
   },
   "outputs": [],
   "source": [
    "# spark = glow.register(spark)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "application/vnd.databricks.v1+cell": {
     "cellMetadata": {},
     "inputWidgets": {},
     "nuid": "1bbdbac5-ab68-4252-8f54-2b59d1eed86e",
     "showTitle": false,
     "title": ""
    }
   },
   "outputs": [],
   "source": [
    "# spark"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "application/vnd.databricks.v1+cell": {
     "cellMetadata": {},
     "inputWidgets": {},
     "nuid": "eb1f1f96-e394-4167-9c51-d8dcca6e9ee5",
     "showTitle": false,
     "title": ""
    },
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# Page Break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "application/vnd.databricks.v1+cell": {
     "cellMetadata": {},
     "inputWidgets": {},
     "nuid": "01e904b3-176a-4c04-8991-bf8c9e313eb9",
     "showTitle": false,
     "title": ""
    },
    "tags": []
   },
   "source": [
    "# Explore bcBio both Somatic and Germline VCFs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "application/vnd.databricks.v1+cell": {
     "cellMetadata": {},
     "inputWidgets": {},
     "nuid": "a980b288-f8d2-45c3-86e5-efd0e0e21183",
     "showTitle": false,
     "title": ""
    },
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "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 #won't work for files in repo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "application/vnd.databricks.v1+cell": {
     "cellMetadata": {},
     "inputWidgets": {},
     "nuid": "7f4d793c-1cbf-4d42-a839-1bf7fc614448",
     "showTitle": false,
     "title": ""
    }
   },
   "outputs": [],
   "source": [
    "%fs\n",
    "ls dbfs:/tmp/dannywong/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "application/vnd.databricks.v1+cell": {
     "cellMetadata": {},
     "inputWidgets": {},
     "nuid": "00a986ed-3d93-4dc9-beb1-fc52f2874d29",
     "showTitle": false,
     "title": ""
    },
    "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\"\n",
    "bcbio_src = \"dbfs:/tmp/dannywong/*.vcf.gz\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "application/vnd.databricks.v1+cell": {
     "cellMetadata": {},
     "inputWidgets": {},
     "nuid": "790effe4-8daa-4f31-8ab1-c44be9267484",
     "showTitle": false,
     "title": ""
    }
   },
   "outputs": [],
   "source": [
    "bcbio_df = spark.read.format(\"vcf\").load(bcbio_src)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "application/vnd.databricks.v1+cell": {
     "cellMetadata": {},
     "inputWidgets": {},
     "nuid": "076fcc40-b741-4f1e-9813-b4e99af6a68a",
     "showTitle": false,
     "title": ""
    }
   },
   "outputs": [],
   "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": 0,
   "metadata": {
    "application/vnd.databricks.v1+cell": {
     "cellMetadata": {},
     "inputWidgets": {},
     "nuid": "890aaa54-2786-4a89-9b61-c3c1a917445c",
     "showTitle": false,
     "title": ""
    }
   },
   "outputs": [],
   "source": [
    "bcbio_df.createOrReplaceTempView(\"vcf_table\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "application/vnd.databricks.v1+cell": {
     "cellMetadata": {},
     "inputWidgets": {},
     "nuid": "874b15da-592c-4d01-865e-82f10d99824b",
     "showTitle": false,
     "title": ""
    }
   },
   "outputs": [],
   "source": [
    "spark.sql(\"describe vcf_table\").show(n=1000, truncate=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "application/vnd.databricks.v1+cell": {
     "cellMetadata": {},
     "inputWidgets": {},
     "nuid": "788b4478-d550-484e-b463-3c4d16f45250",
     "showTitle": false,
     "title": ""
    }
   },
   "outputs": [],
   "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": 0,
   "metadata": {
    "application/vnd.databricks.v1+cell": {
     "cellMetadata": {},
     "inputWidgets": {},
     "nuid": "3abcd40b-8d5f-43c2-8a34-233083bd9d99",
     "showTitle": false,
     "title": ""
    }
   },
   "outputs": [],
   "source": [
    "spark.sql(\"select distinct contigName from vcf_table order by contigName\").show(46)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "application/vnd.databricks.v1+cell": {
     "cellMetadata": {},
     "inputWidgets": {},
     "nuid": "b6e3b254-7d53-49f8-8008-0f75b2f0ff1c",
     "showTitle": false,
     "title": ""
    }
   },
   "outputs": [],
   "source": [
    "spark.sql(\"select contigName, start, end from vcf_table\").show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "application/vnd.databricks.v1+cell": {
     "cellMetadata": {},
     "inputWidgets": {},
     "nuid": "1c2a27a7-3951-4205-849b-ddbacdc3a231",
     "showTitle": false,
     "title": ""
    }
   },
   "outputs": [],
   "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": 0,
   "metadata": {
    "application/vnd.databricks.v1+cell": {
     "cellMetadata": {},
     "inputWidgets": {},
     "nuid": "1f518c80-b6f2-4c54-a0a2-cbbfc4978bab",
     "showTitle": false,
     "title": ""
    }
   },
   "outputs": [],
   "source": [
    "spark.sql(\"select referenceAllele, alternateAlleles, array_size(alternateAlleles) from vcf_table\").show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "application/vnd.databricks.v1+cell": {
     "cellMetadata": {},
     "inputWidgets": {},
     "nuid": "3a4e6d85-547f-46b3-8c13-0f1bf4c52b68",
     "showTitle": false,
     "title": ""
    }
   },
   "outputs": [],
   "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": 0,
   "metadata": {
    "application/vnd.databricks.v1+cell": {
     "cellMetadata": {},
     "inputWidgets": {},
     "nuid": "3b257f71-600f-48f5-ad76-352a1e76051e",
     "showTitle": false,
     "title": ""
    }
   },
   "outputs": [],
   "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": 0,
   "metadata": {
    "application/vnd.databricks.v1+cell": {
     "cellMetadata": {},
     "inputWidgets": {},
     "nuid": "7425d08c-9a1d-4d38-a4d9-46ac70936a98",
     "showTitle": false,
     "title": ""
    }
   },
   "outputs": [],
   "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": 0,
   "metadata": {
    "application/vnd.databricks.v1+cell": {
     "cellMetadata": {},
     "inputWidgets": {},
     "nuid": "880963ac-4e7e-4a8d-ab27-d1508adc4014",
     "showTitle": false,
     "title": ""
    }
   },
   "outputs": [],
   "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": 0,
   "metadata": {
    "application/vnd.databricks.v1+cell": {
     "cellMetadata": {},
     "inputWidgets": {},
     "nuid": "0e80e48b-2f52-420e-a50b-e899bcc0eb66",
     "showTitle": false,
     "title": ""
    },
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# Page Break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "application/vnd.databricks.v1+cell": {
     "cellMetadata": {},
     "inputWidgets": {},
     "nuid": "81a52b7f-eba8-402e-875c-3007ac330ea7",
     "showTitle": false,
     "title": ""
    }
   },
   "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": 0,
   "metadata": {
    "application/vnd.databricks.v1+cell": {
     "cellMetadata": {},
     "inputWidgets": {},
     "nuid": "3aa26c2c-1e3a-4640-8e89-61d55006d869",
     "showTitle": false,
     "title": ""
    },
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# Count total records\n",
    "spark.sql(\"select count(1) as number_of_variants from vcf_table\").show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "application/vnd.databricks.v1+cell": {
     "cellMetadata": {},
     "inputWidgets": {},
     "nuid": "3f052fe1-44d9-4e31-a0cf-cb7d36477adc",
     "showTitle": false,
     "title": ""
    },
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "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": 0,
   "metadata": {
    "application/vnd.databricks.v1+cell": {
     "cellMetadata": {},
     "inputWidgets": {},
     "nuid": "a4beea9f-8bf2-48ac-b759-7ec844175ee3",
     "showTitle": false,
     "title": ""
    },
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "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": {
    "application/vnd.databricks.v1+cell": {
     "cellMetadata": {},
     "inputWidgets": {},
     "nuid": "b94131de-2a51-43f0-817c-46bc981e033f",
     "showTitle": false,
     "title": ""
    }
   },
   "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\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "application/vnd.databricks.v1+cell": {
     "cellMetadata": {},
     "inputWidgets": {},
     "nuid": "a893552f-6608-4fda-aa82-5206e993d554",
     "showTitle": false,
     "title": ""
    },
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# Page Break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "application/vnd.databricks.v1+cell": {
     "cellMetadata": {},
     "inputWidgets": {},
     "nuid": "6715e9f7-da51-440f-918f-2b0e818f8ace",
     "showTitle": false,
     "title": ""
    }
   },
   "source": [
    "# Stop Spark Session"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "application/vnd.databricks.v1+cell": {
     "cellMetadata": {},
     "inputWidgets": {},
     "nuid": "89326173-5dd2-44f8-bd00-749e9e116d38",
     "showTitle": false,
     "title": ""
    },
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "spark.stop()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "application/vnd.databricks.v1+cell": {
     "cellMetadata": {},
     "inputWidgets": {},
     "nuid": "3e17ccb7-f664-431f-b7f2-12b8f466f4f0",
     "showTitle": false,
     "title": ""
    },
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# Continue to next notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "application/vnd.databricks.v1+cell": {
     "cellMetadata": {},
     "inputWidgets": {},
     "nuid": "6639561b-0931-4bdc-818d-7fb205c5ad85",
     "showTitle": false,
     "title": ""
    },
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "application/vnd.databricks.v1+notebook": {
   "dashboards": [],
   "language": "python",
   "notebookMetadata": {
    "mostRecentlyExecutedCommandWithImplicitDF": {
     "commandId": 649139525341007,
     "dataframes": [
      "_sqldf"
     ]
    },
    "pythonIndentUnit": 2
   },
   "notebookName": "03_glow_vcf_explore_merged_databricks",
   "notebookOrigID": 649139525340952,
   "widgets": {}
  },
  "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": 0
}