From dd4d1d7b35cc7b71d96b964e444f6b8ee4bb3a38 Mon Sep 17 00:00:00 2001 From: jsture Date: Fri, 19 Dec 2025 22:12:53 +0100 Subject: [PATCH 1/3] shameful huge commit complete overhaul --- notebooks/WGS/1_QC_WGS.ipynb | 1265 ++++++++++++++++++---------------- 1 file changed, 670 insertions(+), 595 deletions(-) diff --git a/notebooks/WGS/1_QC_WGS.ipynb b/notebooks/WGS/1_QC_WGS.ipynb index 65a880d..55a19c8 100644 --- a/notebooks/WGS/1_QC_WGS.ipynb +++ b/notebooks/WGS/1_QC_WGS.ipynb @@ -1,597 +1,672 @@ { - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# SCRIPT TO PERFORM QUALITY CONTROL ON WHOLE-GENOME SEQUENCING DATA\n", - "\n", - "In order to run, there has to be several files in the project folder:\n", - "- GENCODE GTF: Run Scripts/WGS/01_get_gencode_annotation.sh. Obtain from: https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_46/gencode.v46.annotation.gtf.gz (Check for newer versions).\n", - "\n", - "Once completed, a new Jupyter Notebook should be initialized so we can access this file. Or unmount and mount again the project (https://community.ukbiobank.ac.uk/hc/en-gb/community/posts/16019592366365-It-seems-that-the-recently-dx-uploaded-files-does-not-show-up-on-mnt-project-until-I-re-start-the-whole-Jupyter-Lab-VM)\n", - "\n", - "\n", - "- PVCF BLOCKS: Run Notebooks/WGS/DragenBlockProcessing.ipynb. Obtain from: https://biobank.ndph.ox.ac.uk/ukb/ukb/auxdata/dragen_pvcf_coordinates.zip \n", - "It needs parsing, but in data/misc it is already parsed.\n", - "\n", - "- Samples to remove file: Run Notebooks/WGS/01_QC_Samples.ipynb\n", - "\n", - "#### Initialization \n", - "##### Load packages\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import dxpy\n", - "import pyspark\n", - "\n", - "import hail as hl\n", - "from pathlib import Path\n", - "from datetime import datetime\n", - "\n", - "from src.matrixtables import smart_split_multi_mt" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Constants\n", - "DATABASE = \"matrix_tables\"\n", - "REFERENCE_GENOME = \"GRCh38\"\n", - "PROJ_NAME = \"GIPR\"\n", - "\n", - "Path(\"/tmp\").resolve().mkdir(parents=True, exist_ok=True)\n", - "\n", - "LOG_FILE = (\n", - " Path(\"../hail_logs\", f\"{PROJ_NAME}_{datetime.now().strftime('%H%M')}.log\")\n", - " .resolve()\n", - " .__str__()\n", - ")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Hail and spark configuration" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Spark init\n", - "sc = pyspark.SparkContext()\n", - "spark = pyspark.sql.SparkSession(sc)\n", - "\n", - "# Create database in DNAX\n", - "spark.sql(f\"CREATE DATABASE IF NOT EXISTS {DATABASE} LOCATION 'dnax://'\")\n", - "mt_database = dxpy.find_one_data_object(name=DATABASE)[\"id\"]\n", - "\n", - "# Hail init\n", - "hl.init(sc=sc, default_reference=REFERENCE_GENOME, log=LOG_FILE)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Variables" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# RAP\n", - "VCF_VERSION = \"v1\"\n", - "FIELD_ID = 24310 # DRAGEN population level WGS variants, pVCF format 500k release\n", - "\n", - "# Paths\n", - "BULK_DIR = Path(\"/mnt/project/Bulk\")\n", - "\n", - "# Genes\n", - "GENES = [\"GIPR\"]" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Quality control\n", - "\n", - "#### Gene intervals and blocks " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Get gene intervals\n", - "gene_interval = hl.experimental.get_gene_intervals(\n", - " gene_symbols=GENES,\n", - " reference_genome=\"GRCh38\",\n", - " gtf_file=\"file:///mnt/project/WGS_Javier/WGS_QC/gencode.v46.annotation.gtf\",\n", - ")\n", - "gene_interval" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Get DRAGEN pVCF blocks\n", - "blocks = hl.import_table(\"file:///mnt/project/WGS_Javier/WGS_QC/dragen_pvcf_blocks.tsv\", no_header=False)\n", - "blocks = blocks.annotate(Chromosome=blocks.Chromosome.replace(\"23\", \"X\").replace(\"24\", \"Y\"))\n", - "blocks = blocks.annotate(region=hl.str(\"\").join([hl.str(\"chr\"), blocks.Chromosome]))\n", - "blocks = blocks.annotate(\n", - " interval=hl.locus_interval(\n", - " blocks.region,\n", - " hl.int32(blocks.Starting_Position),\n", - " hl.int32(blocks.Ending_Position),\n", - " reference_genome=\"GRCh38\",\n", - " )\n", - ").key_by(\"interval\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Get blocks for given genes\n", - "gb = blocks.filter(hl.any(lambda inter: blocks.interval.overlaps(inter), gene_interval))\n", - "gb.show()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Import vcf files of specific blocks" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "VCF_DIR = Path(\"DRAGEN WGS/DRAGEN population level WGS variants, pVCF format 500k release\")\n", - "\n", - "vcf_files = [\n", - " f\"file://{BULK_DIR / VCF_DIR}/{chromosome}/ukb{FIELD_ID}_c{chromosome.replace('chr', '')}_b{block}_{VCF_VERSION}.vcf.gz\"\n", - " for block, chromosome in zip(gb.Block.collect(), gb.region.collect())\n", - "]\n", - "\n", - "mt = hl.import_vcf(\n", - " vcf_files,\n", - " drop_samples=False,\n", - " reference_genome=\"GRCh38\",\n", - " array_elements_required=False,\n", - " force_bgz=True,\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Only genes of interest\n", - "mt = hl.filter_intervals(mt, gene_interval)\n", - "print(f\"{mt.count_rows()} variants after gene filtering\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# First checkpoint\n", - "stage = \"FIRST\"\n", - "checkpoint_file = f\"/tmp/{PROJ_NAME}.{stage}.cp.mt\"\n", - "\n", - "mt = mt.checkpoint(checkpoint_file, overwrite=True)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Multi-allele filtering" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Remove variants with 6 or more alleles\n", - "mt = mt.filter_rows(mt.alleles.length() <= 6)\n", - "print(f\"{mt.count_rows()} variants with not more than 6 alleles\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Split multi-allele variants into single ones\n", - "mt = smart_split_multi_mt(mt)\n", - "print(f\"{mt.count_rows()} variants after multi-allele splitting\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Quality control filtering" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "mt = mt.filter_entries(mt.FT == \"PASS\")\n", - "\n", - "# Then, filter variants where there is at least one non-missing genotype\n", - "mt = mt.filter_rows(hl.agg.any(hl.is_defined(mt.GT)))\n", - "print(f\"{mt.count_rows()} variants after entries quality filtering\")\n", - "print(f\"{mt.count_cols()} samples after entries quality filtering\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Compute statistics about the number and fraction of filtered entries.\n", - "mt = hl.MatrixTable.compute_entry_filter_stats(mt, row_field='entry_stats_row', col_field='entry_stats_col')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "row_fraction_threshold = 0.95\n", - "\n", - "# Filter variants where at least 95% of genotypes are unfiltered\n", - "mt = mt.filter_rows(\n", - " (1 - mt.entry_stats_row.fraction_filtered) > row_fraction_threshold\n", - ")\n", - "\n", - "print(f\"{mt.count_rows()} variants after fraction-based filtering\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "col_fraction_threshold = 0.95\n", - "\n", - "# Filter samples where at least 95% of variants are unfiltered\n", - "mt = mt.filter_cols(\n", - " (1 - mt.entry_stats_col.fraction_filtered) > col_fraction_threshold\n", - ")\n", - "\n", - "print(f\"{mt.count_cols()} samples after fraction-based filtering\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Remove samples from 01_QC_Samples.ipynb" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "samples_to_remove = hl.import_table(\"file:///mnt/project/WGS_Javier/Data/Input_regenie/samples_to_remove.tsv\", key=\"eid\")\n", - "\n", - "mt = mt.anti_join_cols(samples_to_remove)\n", - "\n", - "# Filter rows (variants) where any sample information is still present\n", - "mt = mt.filter_rows(hl.agg.any(mt.GT.n_alt_alleles() > 0))\n", - "\n", - "print(f\"Samples remaining after removing samples from QC samples: {mt.count_cols()} \")\n", - "print(f\"Variants remaining after removing samples from QC samples: {mt.count_rows()} \")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Second checkpoint\n", - "stage = \"SECOND\"\n", - "checkpoint_file = f\"/tmp/{PROJ_NAME}.{stage}.cp.mt\"\n", - "\n", - "mt = mt.checkpoint(checkpoint_file, overwrite=True)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Variant Effect Predictor (VEP)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "VEP_JSON = Path(\"GRCh38_VEP.json\").resolve()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "mt = hl.vep(mt, f\"file:{VEP_JSON}\", block_size = 100)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "is_MANE = mt.aggregate_rows(\n", - " hl.agg.all(hl.is_defined(mt.vep.transcript_consequences.mane_select))\n", - ")\n", - "assert is_MANE, \"Selected transcript may not be MANE Select. Check manually.\"\n", - "\n", - "mt = mt.annotate_rows(\n", - " protCons=mt.vep.transcript_consequences.amino_acids[0].split(\"/\")[0]\n", - " + hl.str(mt.vep.transcript_consequences.protein_end[0])\n", - " + mt.vep.transcript_consequences.amino_acids[0].split(\"/\")[-1],\n", - " varid=hl.variant_str(mt.locus, mt.alleles)\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Third checkpoint\n", - "stage = \"THIRD\"\n", - "checkpoint_file = f\"/tmp/{PROJ_NAME}.{stage}.cp.mt\"\n", - "\n", - "mt = mt.checkpoint(checkpoint_file, overwrite=True)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Filtering" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "mt = hl.variant_qc(mt)\n", - "mt = mt.filter_rows(mt.variant_qc.n_non_ref > 0)\n", - "mt = mt.filter_rows(mt.variant_qc.gq_stats.mean >= 20)\n", - "mt = mt.filter_rows(mt.variant_qc.call_rate >= 0.95)\n", - "mt = mt.filter_rows(mt.vep.most_severe_consequence != \"intron_variant\")\n", - "mt = mt.filter_rows(mt.vep.most_severe_consequence != \"downstream_gene_variant\")\n", - "mt = mt.filter_rows(mt.vep.most_severe_consequence != \"upstream_gene_variant\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "print(f\"{mt.count_rows()} variants after quality filtering\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Forth checkpoint\n", - "stage = \"FORTH\"\n", - "checkpoint_file = f\"/tmp/{PROJ_NAME}.{stage}.cp.mt\"\n", - "\n", - "mt = mt.checkpoint(checkpoint_file, overwrite=True)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Formating" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "qt = mt.rows()\n", - "qt = qt.explode(qt.vep.transcript_consequences)\n", - "\n", - "qt = qt.select(\n", - " qt.varid,\n", - " qt.protCons,\n", - " qt.vep.most_severe_consequence,\n", - " qt.vep.transcript_consequences.protein_end,\n", - " qt.vep.transcript_consequences.protein_start,\n", - " qt.vep.transcript_consequences.amino_acids,\n", - " qt.vep.transcript_consequences.gene_id,\n", - " qt.vep.transcript_consequences.transcript_id,\n", - " **qt.variant_qc.flatten(),\n", - ")\n", - "\n", - "qt = qt.annotate(AC=qt.AC[1], AF=qt.AF[1], homozygote_count=qt.homozygote_count[1])\n", - "qt = qt.key_by().drop(\"locus\", \"alleles\")\n", - "\n", - "qt.show(5)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Group by each distinct 'most_severe_consequence' and count the number of occurrences\n", - "consequence_counts = qt.aggregate(\n", - " hl.agg.group_by(qt.most_severe_consequence, hl.agg.count())\n", - ")\n", - "\n", - "print(consequence_counts)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Export " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "qt.export(\"/tmp/variant_qc.tsv\")\n", - "!hadoop fs -getmerge /tmp/variant_qc.tsv ../variant_qc.tsv\n", - "!dx upload ../variant_qc.tsv --path /WGS_Javier/WGS_QC/Output/" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# BGEN file\n", - "BGEN_FILE = \"/tmp/GIPR\"\n", - "GPs = hl.literal([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]])\n", - "\n", - "mt = mt.annotate_entries(GP=GPs[mt.GT.n_alt_alleles()])\n", - "\n", - "hl.export_bgen(\n", - " mt=mt, varid=mt.varid, rsid=mt.varid, gp=mt.GP, output=\"file:\" + BGEN_FILE\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# ANNOTATIONS file\n", - "ANNOTATIONS_FILE = \"/tmp/GIPR.annotations\"\n", - "\n", - "annotations = (\n", - " mt.select_rows(\n", - " varid=mt.varid,\n", - " gene=mt.vep.transcript_consequences.gene_symbol[0],\n", - " annotation=hl.if_else(\n", - " # Check if 'protCons' is missing, if so, use \"most_severe_consequence\"\n", - " hl.is_missing(mt.protCons), \n", - " mt.vep.most_severe_consequence, \n", - " mt.protCons \n", - " )\n", - " )\n", - " .rows()\n", - " .key_by(\"varid\")\n", - " .drop(\"locus\")\n", - " .drop(\"alleles\")\n", - ")\n", - "\n", - "annotations.export(\"file:\" + ANNOTATIONS_FILE, header=False)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# SETLIST file\n", - "SETLIST_FILE = \"/tmp/GIPR.setlist\"\n", - "position = mt.aggregate_rows(hl.agg.min(mt.locus.position))\n", - "names = mt.varid.collect()\n", - "names_str = \",\".join(names)\n", - "\n", - "line = f\"{mt.vep.transcript_consequences.gene_symbol.collect()[0]}\\t{mt.locus.contig.collect()[0]}\\t{position}\\t{names_str}\"\n", - "\n", - "with open(SETLIST_FILE, \"w\") as f:\n", - " f.write(line)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "bgen_file = BGEN_FILE + \".bgen\"\n", - "sample_file = BGEN_FILE + \".sample\"\n", - "\n", - "!dx upload $bgen_file $sample_file $ANNOTATIONS_FILE $SETLIST_FILE --path /WGS_Javier/WGS_QC/Output/" - ] - } - ], - "metadata": { - "language_info": { - "name": "python" - } - }, - "nbformat": 4, - "nbformat_minor": 2 + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# WGS QUALITY CONTROL & PREPARATION PIPELINE\n", + "\n", + "### Overview\n", + "This notebook performs extraction, quality control, and formatting of UK Biobank Whole Genome Sequencing (WGS) data for specific genes of interest. It is optimized for the **DRAGEN 500k Release** (pVCF format).\n", + "\n", + "### VEP is gone\n", + "VEP functionality has been removed from this notebook since it was brittle and heavy. Currently, .annotation and .setlist files have to be built manually until someone steps up and vibe codes an actual notebook based on the exported _variants.tsv\n", + "\n", + "### Workflow Logic\n", + "1. **Target Selection:** Loads pre-computed gene coordinates and identifies the specific VCF data blocks containing these genes.\n", + "2. **Data Extraction:** Imports only the relevant VCF blocks and immediately filters to the gene interval to minimize memory usage.\n", + "3. **Sample Quality Control:** Removes samples identified as \"poor quality\" or \"related\" from the upstream QC notebook (*01_QC_Samples*).\n", + "4. **Variant Quality Control:** * Filters for `FT == PASS` (or missing, as per DRAGEN standards).\n", + " * Removes variants with low call rates (<95%) or no alternate alleles.\n", + "5. **Export:**\n", + " * **Genotypes (.bgen):** Converts Phred-scaled likelihoods (PL) to Genotype Probabilities (GP) and exports to BGEN format for association testing (e.g., Regenie/SAIGE).\n", + " * **Statistics (.tsv):** Exports a clean summary of passing variants (AF, AC, Call Rate) for record-keeping.\n", + "\n", + "---\n", + "\n", + "### Required Input Files\n", + "Before running, ensure the following files are present in the project. If you have just uploaded them, you may need to restart the Jupyter Lab VM to see them.\n", + "\n", + "#### 1. Gene-VCF Overlaps File (`gene_vcf_overlaps.tsv`)\n", + "* **Description:** A mapping file that links Gene Symbols to their genomic coordinates and the specific UKB VCF block(s) that contain them.\n", + "* **Source:** Generated by mapping the GENCODE GTF against the UKB VCF block coordinates.\n", + "* **Columns:** `gene_name`, `chromosome`, `start`, `stop`, `overlapping_vcfs`.\n", + "\n", + "#### 2. Samples to Remove (`samples_to_remove.tsv`)\n", + "* **Description:** A list of Sample IDs (EIDs) to exclude from analysis.\n", + "* **Source:** Output of `Notebooks/WGS/01_QC_Samples.ipynb`.\n", + "* **Content:** Includes withdrawn consent (handled by UKB/DNAx automatically, but good to verify), high-missingness samples, sex discordance, and related individuals (pruned to maximal independent set).\n", + "\n", + "---\n", + "\n", + "#### Initialization \n", + "##### Load packages" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + "
\n", + " \n", + " Loading BokehJS ...\n", + "
\n" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/javascript": "(function(root) {\n function now() {\n return new Date();\n }\n\n const force = true;\n\n if (typeof root._bokeh_onload_callbacks === \"undefined\" || force === true) {\n root._bokeh_onload_callbacks = [];\n root._bokeh_is_loading = undefined;\n }\n\nconst JS_MIME_TYPE = 'application/javascript';\n const HTML_MIME_TYPE = 'text/html';\n const EXEC_MIME_TYPE = 'application/vnd.bokehjs_exec.v0+json';\n const CLASS_NAME = 'output_bokeh rendered_html';\n\n /**\n * Render data to the DOM node\n */\n function render(props, node) {\n const script = document.createElement(\"script\");\n node.appendChild(script);\n }\n\n /**\n * Handle when an output is cleared or removed\n */\n function handleClearOutput(event, handle) {\n const cell = handle.cell;\n\n const id = cell.output_area._bokeh_element_id;\n const server_id = cell.output_area._bokeh_server_id;\n // Clean up Bokeh references\n if (id != null && id in Bokeh.index) {\n Bokeh.index[id].model.document.clear();\n delete Bokeh.index[id];\n }\n\n if (server_id !== undefined) {\n // Clean up Bokeh references\n const cmd_clean = \"from bokeh.io.state import curstate; print(curstate().uuid_to_server['\" + server_id + \"'].get_sessions()[0].document.roots[0]._id)\";\n cell.notebook.kernel.execute(cmd_clean, {\n iopub: {\n output: function(msg) {\n const id = msg.content.text.trim();\n if (id in Bokeh.index) {\n Bokeh.index[id].model.document.clear();\n delete Bokeh.index[id];\n }\n }\n }\n });\n // Destroy server and session\n const cmd_destroy = \"import bokeh.io.notebook as ion; ion.destroy_server('\" + server_id + \"')\";\n cell.notebook.kernel.execute(cmd_destroy);\n }\n }\n\n /**\n * Handle when a new output is added\n */\n function handleAddOutput(event, handle) {\n const output_area = handle.output_area;\n const output = handle.output;\n\n // limit handleAddOutput to display_data with EXEC_MIME_TYPE content only\n if ((output.output_type != \"display_data\") || (!Object.prototype.hasOwnProperty.call(output.data, EXEC_MIME_TYPE))) {\n return\n }\n\n const toinsert = output_area.element.find(\".\" + CLASS_NAME.split(' ')[0]);\n\n if (output.metadata[EXEC_MIME_TYPE][\"id\"] !== undefined) {\n toinsert[toinsert.length - 1].firstChild.textContent = output.data[JS_MIME_TYPE];\n // store reference to embed id on output_area\n output_area._bokeh_element_id = output.metadata[EXEC_MIME_TYPE][\"id\"];\n }\n if (output.metadata[EXEC_MIME_TYPE][\"server_id\"] !== undefined) {\n const bk_div = document.createElement(\"div\");\n bk_div.innerHTML = output.data[HTML_MIME_TYPE];\n const script_attrs = bk_div.children[0].attributes;\n for (let i = 0; i < script_attrs.length; i++) {\n toinsert[toinsert.length - 1].firstChild.setAttribute(script_attrs[i].name, script_attrs[i].value);\n toinsert[toinsert.length - 1].firstChild.textContent = bk_div.children[0].textContent\n }\n // store reference to server id on output_area\n output_area._bokeh_server_id = output.metadata[EXEC_MIME_TYPE][\"server_id\"];\n }\n }\n\n function register_renderer(events, OutputArea) {\n\n function append_mime(data, metadata, element) {\n // create a DOM node to render to\n const toinsert = this.create_output_subarea(\n metadata,\n CLASS_NAME,\n EXEC_MIME_TYPE\n );\n this.keyboard_manager.register_events(toinsert);\n // Render to node\n const props = {data: data, metadata: metadata[EXEC_MIME_TYPE]};\n render(props, toinsert[toinsert.length - 1]);\n element.append(toinsert);\n return toinsert\n }\n\n /* Handle when an output is cleared or removed */\n events.on('clear_output.CodeCell', handleClearOutput);\n events.on('delete.Cell', handleClearOutput);\n\n /* Handle when a new output is added */\n events.on('output_added.OutputArea', handleAddOutput);\n\n /**\n * Register the mime type and append_mime function with output_area\n */\n OutputArea.prototype.register_mime_type(EXEC_MIME_TYPE, append_mime, {\n /* Is output safe? */\n safe: true,\n /* Index of renderer in `output_area.display_order` */\n index: 0\n });\n }\n\n // register the mime type if in Jupyter Notebook environment and previously unregistered\n if (root.Jupyter !== undefined) {\n const events = require('base/js/events');\n const OutputArea = require('notebook/js/outputarea').OutputArea;\n\n if (OutputArea.prototype.mime_types().indexOf(EXEC_MIME_TYPE) == -1) {\n register_renderer(events, OutputArea);\n }\n }\n if (typeof (root._bokeh_timeout) === \"undefined\" || force === true) {\n root._bokeh_timeout = Date.now() + 5000;\n root._bokeh_failed_load = false;\n }\n\n const NB_LOAD_WARNING = {'data': {'text/html':\n \"
\\n\"+\n \"

\\n\"+\n \"BokehJS does not appear to have successfully loaded. If loading BokehJS from CDN, this \\n\"+\n \"may be due to a slow or bad network connection. Possible fixes:\\n\"+\n \"

\\n\"+\n \"\\n\"+\n \"\\n\"+\n \"from bokeh.resources import INLINE\\n\"+\n \"output_notebook(resources=INLINE)\\n\"+\n \"\\n\"+\n \"
\"}};\n\n function display_loaded() {\n const el = document.getElementById(\"p1001\");\n if (el != null) {\n el.textContent = \"BokehJS is loading...\";\n }\n if (root.Bokeh !== undefined) {\n if (el != null) {\n el.textContent = \"BokehJS \" + root.Bokeh.version + \" successfully loaded.\";\n }\n } else if (Date.now() < root._bokeh_timeout) {\n setTimeout(display_loaded, 100)\n }\n }\n\n function run_callbacks() {\n try {\n root._bokeh_onload_callbacks.forEach(function(callback) {\n if (callback != null)\n callback();\n });\n } finally {\n delete root._bokeh_onload_callbacks\n }\n console.debug(\"Bokeh: all callbacks have finished\");\n }\n\n function load_libs(css_urls, js_urls, callback) {\n if (css_urls == null) css_urls = [];\n if (js_urls == null) js_urls = [];\n\n root._bokeh_onload_callbacks.push(callback);\n if (root._bokeh_is_loading > 0) {\n console.debug(\"Bokeh: BokehJS is being loaded, scheduling callback at\", now());\n return null;\n }\n if (js_urls == null || js_urls.length === 0) {\n run_callbacks();\n return null;\n }\n console.debug(\"Bokeh: BokehJS not loaded, scheduling load and callback at\", now());\n root._bokeh_is_loading = css_urls.length + js_urls.length;\n\n function on_load() {\n root._bokeh_is_loading--;\n if (root._bokeh_is_loading === 0) {\n console.debug(\"Bokeh: all BokehJS libraries/stylesheets loaded\");\n run_callbacks()\n }\n }\n\n function on_error(url) {\n console.error(\"failed to load \" + url);\n }\n\n for (let i = 0; i < css_urls.length; i++) {\n const url = css_urls[i];\n const element = document.createElement(\"link\");\n element.onload = on_load;\n element.onerror = on_error.bind(null, url);\n element.rel = \"stylesheet\";\n element.type = \"text/css\";\n element.href = url;\n console.debug(\"Bokeh: injecting link tag for BokehJS stylesheet: \", url);\n document.body.appendChild(element);\n }\n\n for (let i = 0; i < js_urls.length; i++) {\n const url = js_urls[i];\n const element = document.createElement('script');\n element.onload = on_load;\n element.onerror = on_error.bind(null, url);\n element.async = false;\n element.src = url;\n console.debug(\"Bokeh: injecting script tag for BokehJS library: \", url);\n document.head.appendChild(element);\n }\n };\n\n function inject_raw_css(css) {\n const element = document.createElement(\"style\");\n element.appendChild(document.createTextNode(css));\n document.body.appendChild(element);\n }\n\n const js_urls = [\"https://cdn.bokeh.org/bokeh/release/bokeh-3.0.3.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-gl-3.0.3.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-widgets-3.0.3.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-tables-3.0.3.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-mathjax-3.0.3.min.js\"];\n const css_urls = [];\n\n const inline_js = [ function(Bokeh) {\n Bokeh.set_log_level(\"info\");\n },\nfunction(Bokeh) {\n }\n ];\n\n function run_inline_js() {\n if (root.Bokeh !== undefined || force === true) {\n for (let i = 0; i < inline_js.length; i++) {\n inline_js[i].call(root, root.Bokeh);\n }\nif (force === true) {\n display_loaded();\n }} else if (Date.now() < root._bokeh_timeout) {\n setTimeout(run_inline_js, 100);\n } else if (!root._bokeh_failed_load) {\n console.log(\"Bokeh: BokehJS failed to load within specified timeout.\");\n root._bokeh_failed_load = true;\n } else if (force !== true) {\n const cell = $(document.getElementById(\"p1001\")).parents('.cell').data().cell;\n cell.output_area.append_execute_result(NB_LOAD_WARNING)\n }\n }\n\n if (root._bokeh_is_loading === 0) {\n console.debug(\"Bokeh: BokehJS loaded, going straight to plotting\");\n run_inline_js();\n } else {\n load_libs(css_urls, js_urls, function() {\n console.debug(\"Bokeh: BokehJS plotting callback run at\", now());\n run_inline_js();\n });\n }\n}(window));", + "application/vnd.bokehjs_load.v0+json": "" + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "import dxpy\n", + "import pyspark\n", + "\n", + "import hail as hl\n", + "from pathlib import Path\n", + "from datetime import datetime" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "# Constants\n", + "DATABASE = \"matrix_tables\"\n", + "REFERENCE_GENOME = \"GRCh38\"\n", + "PROJ_NAME = \"GIPR\"\n", + "\n", + "Path(\"/tmp\").resolve().mkdir(parents=True, exist_ok=True)\n", + "\n", + "LOG_FILE = (\n", + " Path(\"../hail_logs\", f\"{PROJ_NAME}_{datetime.now().strftime('%H%M')}.log\")\n", + " .resolve()\n", + " .__str__()\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Hail and spark configuration" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "pip-installed Hail requires additional configuration options in Spark referring\n", + " to the path to the Hail Python module directory HAIL_DIR,\n", + " e.g. /path/to/python/site-packages/hail:\n", + " spark.jars=HAIL_DIR/backend/hail-all-spark.jar\n", + " spark.driver.extraClassPath=HAIL_DIR/backend/hail-all-spark.jar\n", + " spark.executor.extraClassPath=./hail-all-spark.jarRunning on Apache Spark version 3.5.2\n", + "SparkUI available at http://ip-10-60-122-101.eu-west-2.compute.internal:8081\n", + "Welcome to\n", + " __ __ <>__\n", + " / /_/ /__ __/ /\n", + " / __ / _ `/ / /\n", + " /_/ /_/\\_,_/_/_/ version 0.2.132-678e1f52b999\n", + "LOGGING: writing to /opt/hail_logs/GIPR_1917.log\n", + "2025-12-19 19:18:23.509 Hail: INFO: Reading table to impute column types\n", + "2025-12-19 19:18:24.694 Hail: INFO: Finished type imputation\n", + " Loading field 'gene_id' as type str (imputed)\n", + " Loading field 'gene_name' as type str (imputed)\n", + " Loading field 'chromosome' as type str (imputed)\n", + " Loading field 'start' as type int32 (imputed)\n", + " Loading field 'stop' as type int32 (imputed)\n", + " Loading field 'vcf_count' as type int32 (imputed)\n", + " Loading field 'overlapping_vcfs' as type str (imputed)\n", + " Loading field 'vcf_start' as type str (imputed)\n", + " Loading field 'vcf_stop' as type str (imputed)\n", + "2025-12-19 19:19:34.255 Hail: INFO: scanning VCF for sortedness...\n", + "2025-12-19 19:20:03.851 Hail: INFO: Coerced prefix-sorted VCF, requiring additional sorting within data partitions on each query.\n", + "2025-12-19 19:36:59.603 Hail: INFO: wrote matrix table with 16066 rows and 4810 columns in 4 partitions to /tmp/GIPR.FIRST.cp.mt\n", + "2025-12-19 19:55:54.179 Hail: INFO: wrote matrix table with 15667 rows and 490541 columns in 4 partitions to /tmp/GIPR.FIRST.cp.mt\n", + "2025-12-19 19:55:55.453 Hail: INFO: Reading table without type imputation\n", + " Loading field 'eid' as type str (not specified)\n", + "2025-12-19 20:13:48.862 Hail: INFO: wrote matrix table with 13366 rows and 442704 columns in 4 partitions to /tmp/GIPR.SECOND.cp.mt\n", + "2025-12-19 20:19:54.238 Hail: INFO: merging 5 files totalling 1.2M...\n", + "2025-12-19 20:26:56.454 Hail: INFO: merging 5 files totalling 1.2M...\n", + "2025-12-19 20:28:06.649 Hail: INFO: merging 5 files totalling 1.2M...\n", + "2025-12-19 20:28:06.697 Hail: INFO: while writing:\n", + " file:/tmp/GIPR_variants.tsv\n", + " merge time: 47.916ms\n", + "2025-12-19 20:43:56.335 Hail: INFO: merging 5 files totalling 1.3M...\n", + "2025-12-19 20:43:56.353 Hail: INFO: while writing:\n", + " file:/tmp/GIPR_variants.tsv\n", + " merge time: 17.422ms\n", + "2025-12-19 20:52:36.890 Hail: INFO: merging 5 files totalling 1.3M...\n", + "2025-12-19 20:52:36.906 Hail: INFO: while writing:\n", + " file:///tmp/GIPR_variants.tsv\n", + " merge time: 16.407ms\n", + "2025-12-19 21:04:32.893 Hail: INFO: merging 5 files totalling 1.3M...\n", + "2025-12-19 21:04:32.911 Hail: INFO: while writing:\n", + " file:///tmp/GIPR_variants.tsv\n", + " merge time: 17.302ms\n" + ] + } + ], + "source": [ + "# Spark init\n", + "sc = pyspark.SparkContext()\n", + "spark = pyspark.sql.SparkSession(sc)\n", + "\n", + "# Create database in DNAX\n", + "spark.sql(f\"CREATE DATABASE IF NOT EXISTS {DATABASE} LOCATION 'dnax://'\")\n", + "mt_database = dxpy.find_one_data_object(name=DATABASE)[\"id\"]\n", + "\n", + "# Hail init\n", + "hl.init(sc=sc, log=LOG_FILE)\n", + "hl.default_reference(REFERENCE_GENOME)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Variables" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "# RAP\n", + "VCF_VERSION = \"v1\"\n", + "FIELD_ID = (\n", + " 24311 # ML-corrected DRAGEN population level WGS variants, pVCF format 500k release\n", + ")\n", + "VCF_DIR = Path(\n", + " \"DRAGEN WGS/ML-corrected DRAGEN population level WGS variants, pVCF format 500k release\"\n", + ")\n", + "\n", + "# Paths\n", + "BULK_DIR = Path(\"/mnt/project/Bulk\")\n", + "\n", + "# Genes\n", + "GENES = [\"TAS1R2\", \"GIPR\"]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Load\n", + "\n", + "#### Gene intervals and blocks\n", + "We use the lookup table rather than build block list here" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Loaded coordinates for 2 genes: ['TAS1R2', 'GIPR']\n", + "Identified 4 VCF files containing the genes.\n", + "Sample path: file:///mnt/project/Bulk/DRAGEN WGS/ML-corrected DRAGEN population level WGS variants, pVCF format 500k release/chr1/ukb24311_c1_b941_v1.vcf.gz\n", + "Identified 4 VCF files containing the genes.\n", + "file:///mnt/project/Bulk/DRAGEN WGS/ML-corrected DRAGEN population level WGS variants, pVCF format 500k release/chr1/ukb24311_c1_b941_v1.vcf.gz\n" + ] + } + ], + "source": [ + "# --- REPLACEMENT FOR GENE INTERVALS & BLOCKS ---\n", + "# Load the pre-computed overlap file\n", + "OVERLAP_TSV = \"/mnt/project/WGS_Javier/WGS_QC/gene_vcf_overlaps.tsv\" # Update path\n", + "\n", + "overlaps_ht = hl.import_table(f\"file://{OVERLAP_TSV}\", impute=True) # \u2713\n", + "\n", + "# 1. Filter for the genes of interest\n", + "overlaps_ht = overlaps_ht.filter(hl.literal(GENES).contains(overlaps_ht.gene_name))\n", + "\n", + "# 2. Extract Gene Intervals\n", + "# We use hl.parse_locus_interval to handle \"chr:start-stop\" string format automatically\n", + "intervals_data = overlaps_ht.select(\"chromosome\", \"start\", \"stop\").collect()\n", + "\n", + "gene_intervals = [\n", + " hl.parse_locus_interval(\n", + " f\"{row.chromosome}:{row.start}-{row.stop}\", reference_genome=REFERENCE_GENOME\n", + " )\n", + " for row in intervals_data\n", + "]\n", + "\n", + "print(f\"Loaded coordinates for {len(gene_intervals)} genes: {GENES}\")\n", + "\n", + "# 3. Construct VCF Paths\n", + "# We iterate through the rows to get the correct chromosome folder for each file.\n", + "# Path format: {BULK_DIR}/{VCF_DIR}/{chromosome}/{filename}\n", + "vcf_data = overlaps_ht.select(\"chromosome\", \"overlapping_vcfs\").collect()\n", + "vcf_paths = set()\n", + "\n", + "# Ensure VCF_DIR is a Path object (if defined as string in constants)\n", + "VCF_DIR_PATH = BULK_DIR / Path(VCF_DIR)\n", + "\n", + "for row in vcf_data:\n", + " if not row.overlapping_vcfs:\n", + " continue\n", + "\n", + " # Split \"file1.vcf,file2.vcf\" -> [\"file1.vcf\", \"file2.vcf\"]\n", + " files = row.overlapping_vcfs.split(\",\")\n", + "\n", + " for f in files:\n", + " f = f.strip()\n", + " if not f:\n", + " continue\n", + "\n", + " # Use the 'chromosome' column from the TSV (e.g., \"chr19\") for the folder\n", + " full_path = f\"file://{VCF_DIR_PATH}/{row.chromosome}/{f}\"\n", + " vcf_paths.add(full_path)\n", + "\n", + "vcf_files = list(vcf_paths)\n", + "print(f\"Identified {len(vcf_files)} VCF files containing the genes.\")\n", + "print(\"Sample path:\", vcf_files[0] if vcf_files else \"None\")\n", + "print(f\"Identified {len(vcf_files)} VCF files containing the genes.\")" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "----------------------------------------\n", + "Global fields:\n", + " None\n", + "----------------------------------------\n", + "Column fields:\n", + " 's': str\n", + "----------------------------------------\n", + "Row fields:\n", + " 'locus': locus\n", + " 'alleles': array\n", + " 'rsid': str\n", + " 'qual': float64\n", + " 'filters': set\n", + " 'info': struct {\n", + " AC: array, \n", + " AN: int32, \n", + " NS: int32, \n", + " NS_GT: int32, \n", + " NS_NOGT: int32, \n", + " NS_NODATA: int32, \n", + " IC: float64, \n", + " HWE: array, \n", + " ExcHet: array, \n", + " HWE_CHISQ: float64, \n", + " AF: array\n", + " }\n", + "----------------------------------------\n", + "Entry fields:\n", + " 'GT': call\n", + " 'GQ': int32\n", + " 'LAD': array\n", + " 'FT': str\n", + " 'LPL': array\n", + " 'LAA': array\n", + " 'LAF': array\n", + " 'QL': float64\n", + "----------------------------------------\n", + "Column key: ['s']\n", + "Row key: ['locus', 'alleles']\n", + "----------------------------------------\n" + ] + } + ], + "source": [ + "# --- IMPORT VCFS ---\n", + "mt = hl.import_vcf(\n", + " vcf_files,\n", + " drop_samples=False,\n", + " reference_genome=REFERENCE_GENOME,\n", + " array_elements_required=True,\n", + " force_bgz=True,\n", + ")\n", + "\n", + "# TODO: Make nice\n", + "# Also, watch out\n", + "downsample = False\n", + "if downsample:\n", + " mt = mt.sample_cols(0.01)\n", + "\n", + "# Filter MT to exact gene boundaries (removes flanking regions in the 20kb block)\n", + "mt = hl.filter_intervals(mt, gene_intervals)\n", + "\n", + "# Godlike filtering\n", + "mt = mt.filter_rows(hl.len(mt.filters) == 0)\n", + "\n", + "mt.describe()\n", + "\n", + "# Pre-QC count: 16066 variants, 4810 samples" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Checkpoint MT\n", + "Checkpoint after coarse filtering" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "----------------------------------------\n", + "Global fields:\n", + " None\n", + "----------------------------------------\n", + "Column fields:\n", + " 's': str\n", + "----------------------------------------\n", + "Row fields:\n", + " 'locus': locus\n", + " 'alleles': array\n", + " 'rsid': str\n", + " 'qual': float64\n", + " 'filters': set\n", + " 'info': struct {\n", + " AC: array, \n", + " AN: int32, \n", + " NS: int32, \n", + " NS_GT: int32, \n", + " NS_NOGT: int32, \n", + " NS_NODATA: int32, \n", + " IC: float64, \n", + " HWE: array, \n", + " ExcHet: array, \n", + " HWE_CHISQ: float64, \n", + " AF: array\n", + " }\n", + "----------------------------------------\n", + "Entry fields:\n", + " 'GT': call\n", + " 'GQ': int32\n", + " 'LAD': array\n", + " 'FT': str\n", + " 'LPL': array\n", + " 'LAA': array\n", + " 'LAF': array\n", + " 'QL': float64\n", + "----------------------------------------\n", + "Column key: ['s']\n", + "Row key: ['locus', 'alleles']\n", + "----------------------------------------\n", + "Pre-QC count: 15667 variants, 490541 samples\n" + ] + } + ], + "source": [ + "# First checkpoint\n", + "\n", + "stage = \"FIRST\"\n", + "checkpoint_file = f\"/tmp/{PROJ_NAME}.{stage}.cp.mt\"\n", + "\n", + "mt = mt.checkpoint(checkpoint_file, overwrite=True)\n", + "# mt = hl.read_matrix_table(checkpoint_file)\n", + "\n", + "mt.describe()\n", + "print(f\"Pre-QC count: {mt.count_rows()} variants, {mt.count_cols()} samples\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Quality control filtering\n", + "Remove samples from 01_QC_Samples\n", + "\n", + "Filter to FT (empty? Documentation says PASS? Who knows?)\n", + "\n", + "A lot of the entry fields are populated, so check if ever adding more filters. GQ and LPL are empty for example. \n", + "\n", + "GT is as usual the one we want. \n", + "\n", + "variant_qc is run last to save a little compute" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Running Sample QC...\n", + "Running Variant QC...\n" + ] + } + ], + "source": [ + "# --- QC & FILTERING ---\n", + "\n", + "\n", + "# This should match file from 01_QC_Samples.ipynb\n", + "\n", + "SAMPLES_REMOVE_PATH = \"/mnt/project/TASR/Phenotypes/QC/samples_to_remove.tsv\"\n", + "\n", + "if Path(\n", + " SAMPLES_REMOVE_PATH.replace(\"file://\", \"\")\n", + ").exists() or SAMPLES_REMOVE_PATH.startswith(\"dnax://\"):\n", + " samples_to_remove = hl.import_table(f\"file://{SAMPLES_REMOVE_PATH}\", key=\"eid\")\n", + " mt = mt.anti_join_cols(samples_to_remove)\n", + "else:\n", + " print(\n", + " f\"WARNING: QC file not found at {SAMPLES_REMOVE_PATH}. Skipping sample removal.\"\n", + " )\n", + "\n", + "# 2. Filter Entries (Genotype Level)\n", + "# ML-corrected data is cleaner, but we still enforce PASS and basic depth/quality if available\n", + "# TODO: Find some actual documentation here\n", + "\n", + "\n", + "mt = mt.filter_entries(hl.is_missing(mt.FT) | (mt.FT == \"PASS\"))\n", + "\n", + "\n", + "# 4. Filter Samples (Columns)\n", + "# Remove samples with low call rate (< 95%)\n", + "mt = hl.sample_qc(mt)\n", + "mt = mt.filter_cols(mt.sample_qc.call_rate > 0.95)\n", + "\n", + "\n", + "# 3. Run Variant QC\n", + "mt = hl.variant_qc(mt) # needs this for later too\n", + "mt = mt.filter_rows(\n", + " (mt.variant_qc.n_non_ref > 0) # Must have at least one Alt allele\n", + " & (mt.variant_qc.call_rate > 0.95) # High call rate\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Post-QC count: 13366 variants, 442704 samples\n" + ] + } + ], + "source": [ + "# Second checkpoint\n", + "# Hail still likes checkpointing\n", + "\n", + "stage = \"SECOND\"\n", + "checkpoint_file = f\"/tmp/{PROJ_NAME}.{stage}.cp.mt\"\n", + "\n", + "mt = mt.checkpoint(checkpoint_file, overwrite=True)\n", + "print(f\"Post-QC count: {mt.count_rows()} variants, {mt.count_cols()} samples\")\n", + "mt.describe()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Exports\n", + "\n", + "Again, Regenie checks bgen GP but wants GT\n", + "\n", + "Removed all VEP logic from here. " + ] + }, + { + "cell_type": "code", + "execution_count": 36, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Exporting BGEN to /tmp/GIPR.bgen ... (takes a little while)\n", + "Exporting Stats to /tmp/GIPR_variants.tsv ...\n" + ] + } + ], + "source": [ + "# --- EXPORTS (Cleaned Up) ---\n", + "\n", + "# Define Paths\n", + "LOCAL_TMP = \"/tmp\" # Local container storage\n", + "HDFS_TMP = \"file:///tmp\" # Explicitly pointing to local FS for Hail\n", + "DX_OUTPUT_DIR = \"/TASR/Genotypes\"\n", + "BASE_NAME = f\"{PROJ_NAME}\"\n", + "\n", + "# 1. Prepare Annotations (GP & VarID)\n", + "# Fallback to Hard Calls (since we know PL/LPL are likely missing/empty)\n", + "GPs = hl.literal([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]])\n", + "mt = mt.annotate_entries(GP=GPs[mt.GT.n_alt_alleles()])\n", + "\n", + "mt = mt.annotate_rows(\n", + " varid=hl.delimit(\n", + " [mt.locus.contig, hl.str(mt.locus.position), mt.alleles[0], mt.alleles[1]], \":\"\n", + " )\n", + ")\n", + "\n", + "# 2. Export BGEN (To Local /tmp via \"file://\")\n", + "# \"parallel=None\" ensures we get a single .bgen file, not a directory of shards.\n", + "# We use \"file://\" prefix to force Hail to write to the local node's /tmp,\n", + "# making it instantly accessible for dx upload without hadoop fs -get.\n", + "print(f\"Exporting BGEN to {LOCAL_TMP}/{BASE_NAME}.bgen ... (takes a little while)\")\n", + "\n", + "hl.export_bgen(\n", + " mt=mt,\n", + " output=f\"{HDFS_TMP}/{BASE_NAME}\", # file:///tmp/GIPR\n", + " gp=mt.GP,\n", + " varid=mt.varid,\n", + " rsid=mt.varid,\n", + " parallel=None,\n", + ")\n", + "\n", + "# 3. Export Variant Stats TSV\n", + "rows_ht = mt.rows()\n", + "\n", + "export_ht = rows_ht.select(\n", + " chromosome=rows_ht.locus.contig,\n", + " position=rows_ht.locus.position,\n", + " ref=rows_ht.alleles[0],\n", + " alt=rows_ht.alleles[1],\n", + " rsid=rows_ht.rsid,\n", + " call_rate=rows_ht.variant_qc.call_rate,\n", + " AC=rows_ht.variant_qc.AC[1],\n", + " AF=rows_ht.variant_qc.AF[1],\n", + " AN=rows_ht.variant_qc.AN,\n", + " N_HOM=rows_ht.variant_qc.homozygote_count[1],\n", + " N_HET=rows_ht.variant_qc.n_het,\n", + " N_NC=rows_ht.variant_qc.n_not_called,\n", + ")\n", + "\n", + "tsv_name = f\"{BASE_NAME}_variants.tsv\"\n", + "print(f\"Exporting Stats to {LOCAL_TMP}/{tsv_name} ...\")\n", + "export_ht.export(f\"{HDFS_TMP}/{tsv_name}\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# We now upload directly from /tmp where Hail wrote the files\n", + "!dx upload {LOCAL_TMP}/{BASE_NAME}.bgen {LOCAL_TMP}/{BASE_NAME}.sample {LOCAL_TMP}/{tsv_name} --path {DX_OUTPUT_DIR}" + ] + } + ], + "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.12.9" + } + }, + "nbformat": 4, + "nbformat_minor": 4 } From 5861bb8033c015e3acdce04533fa0e7ca0724078 Mon Sep 17 00:00:00 2001 From: jsture Date: Fri, 19 Dec 2025 23:21:24 +0100 Subject: [PATCH 2/3] refactor for speedup --- notebooks/WGS/1_QC_WGS.ipynb | 88 ++++++++++++++++++++---------------- 1 file changed, 49 insertions(+), 39 deletions(-) diff --git a/notebooks/WGS/1_QC_WGS.ipynb b/notebooks/WGS/1_QC_WGS.ipynb index 55a19c8..82c2f89 100644 --- a/notebooks/WGS/1_QC_WGS.ipynb +++ b/notebooks/WGS/1_QC_WGS.ipynb @@ -197,7 +197,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -209,12 +209,17 @@ "VCF_DIR = Path(\n", " \"DRAGEN WGS/ML-corrected DRAGEN population level WGS variants, pVCF format 500k release\"\n", ")\n", + "SAMPLES_REMOVE_PATH = \"/mnt/project/TASR/Phenotypes/QC/samples_to_remove.tsv\"\n", + "\n", "\n", "# Paths\n", "BULK_DIR = Path(\"/mnt/project/Bulk\")\n", "\n", "# Genes\n", - "GENES = [\"TAS1R2\", \"GIPR\"]" + "GENES = [\"TAS1R2\", \"GIPR\"]\n", + "\n", + "# Downsample for testing\n", + "DOWNSAMPLE_FRACTION: float | None = None # Set to 0.01 for 1% sample, 0.1 for 10%, etc." ] }, { @@ -302,7 +307,24 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "if Path(\n", + " SAMPLES_REMOVE_PATH.replace(\"file://\", \"\")\n", + ").exists() or SAMPLES_REMOVE_PATH.startswith(\"dnax://\"):\n", + " samples_to_remove = hl.import_table(f\"file://{SAMPLES_REMOVE_PATH}\", key=\"eid\")\n", + "else:\n", + " samples_to_remove = None\n", + " print(\n", + " f\"WARNING: QC file not found at {SAMPLES_REMOVE_PATH}. Skipping sample removal.\"\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, "metadata": { "tags": [] }, @@ -364,11 +386,15 @@ " force_bgz=True,\n", ")\n", "\n", - "# TODO: Make nice\n", + "# Avoid write when checkpoint\n", + "if samples_to_remove is not None:\n", + " mt = mt.anti_join_cols(samples_to_remove)\n", + "\n", + "\n", "# Also, watch out\n", - "downsample = False\n", - "if downsample:\n", - " mt = mt.sample_cols(0.01)\n", + "if DOWNSAMPLE_FRACTION is not None:\n", + " print(f\"\u26a0\ufe0f DOWNSAMPLING to {DOWNSAMPLE_FRACTION*100}% of samples for testing\")\n", + " mt = mt.sample_cols(p=DOWNSAMPLE_FRACTION, seed=42)\n", "\n", "# Filter MT to exact gene boundaries (removes flanking regions in the 20kb block)\n", "mt = hl.filter_intervals(mt, gene_intervals)\n", @@ -376,9 +402,15 @@ "# Godlike filtering\n", "mt = mt.filter_rows(hl.len(mt.filters) == 0)\n", "\n", - "mt.describe()\n", + "# 2. Filter Entries (Genotype Level)\n", + "# ML-corrected data is cleaner, but we still enforce PASS and basic depth/quality if available\n", + "# TODO: Find some actual documentation here\n", + "# Minimal load for the |\n", + "\n", + "mt = mt.filter_entries(hl.is_missing(mt.FT) | (mt.FT == \"PASS\"))\n", "\n", - "# Pre-QC count: 16066 variants, 4810 samples" + "\n", + "mt.describe()" ] }, { @@ -473,7 +505,7 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": null, "metadata": {}, "outputs": [ { @@ -488,29 +520,6 @@ "source": [ "# --- QC & FILTERING ---\n", "\n", - "\n", - "# This should match file from 01_QC_Samples.ipynb\n", - "\n", - "SAMPLES_REMOVE_PATH = \"/mnt/project/TASR/Phenotypes/QC/samples_to_remove.tsv\"\n", - "\n", - "if Path(\n", - " SAMPLES_REMOVE_PATH.replace(\"file://\", \"\")\n", - ").exists() or SAMPLES_REMOVE_PATH.startswith(\"dnax://\"):\n", - " samples_to_remove = hl.import_table(f\"file://{SAMPLES_REMOVE_PATH}\", key=\"eid\")\n", - " mt = mt.anti_join_cols(samples_to_remove)\n", - "else:\n", - " print(\n", - " f\"WARNING: QC file not found at {SAMPLES_REMOVE_PATH}. Skipping sample removal.\"\n", - " )\n", - "\n", - "# 2. Filter Entries (Genotype Level)\n", - "# ML-corrected data is cleaner, but we still enforce PASS and basic depth/quality if available\n", - "# TODO: Find some actual documentation here\n", - "\n", - "\n", - "mt = mt.filter_entries(hl.is_missing(mt.FT) | (mt.FT == \"PASS\"))\n", - "\n", - "\n", "# 4. Filter Samples (Columns)\n", "# Remove samples with low call rate (< 95%)\n", "mt = hl.sample_qc(mt)\n", @@ -527,7 +536,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": null, "metadata": {}, "outputs": [ { @@ -541,13 +550,14 @@ "source": [ "# Second checkpoint\n", "# Hail still likes checkpointing\n", + "# TODO: consider re-instating but likely not needed here\n", "\n", - "stage = \"SECOND\"\n", - "checkpoint_file = f\"/tmp/{PROJ_NAME}.{stage}.cp.mt\"\n", + "# stage = \"SECOND\"\n", + "# checkpoint_file = f\"/tmp/{PROJ_NAME}.{stage}.cp.mt\"\n", "\n", - "mt = mt.checkpoint(checkpoint_file, overwrite=True)\n", - "print(f\"Post-QC count: {mt.count_rows()} variants, {mt.count_cols()} samples\")\n", - "mt.describe()" + "# mt = mt.checkpoint(checkpoint_file, overwrite=True)\n", + "# print(f\"Post-QC count: {mt.count_rows()} variants, {mt.count_cols()} samples\")\n", + "# mt.describe()" ] }, { From 712ffd2e89e3f7bfdb1c8df78f468d616c303677 Mon Sep 17 00:00:00 2001 From: jsture Date: Sat, 20 Dec 2025 11:26:02 +0100 Subject: [PATCH 3/3] moved setlist functionality, deprecated old WES functions --- src/matrixtables.py | 83 ++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 82 insertions(+), 1 deletion(-) diff --git a/src/matrixtables.py b/src/matrixtables.py index afe16bf..9ad6ead 100644 --- a/src/matrixtables.py +++ b/src/matrixtables.py @@ -4,6 +4,8 @@ import hail as hl import pandas as pd +import warnings + def import_mt( genes: List[str], @@ -13,6 +15,12 @@ def import_mt( field_id: int = 23157, ) -> hl.matrixtable.MatrixTable: """Maps a (list of) genes to their corresponding VCF file and region before import as Hail MatrixTable + + .. deprecated:: + This function is hard deprecated and will be removed in a future version. + hl.import_gvcfs() is no longer available in recent Hail versions. + See WGS 01_QC notebook implementation for alternatives. + Parameters ---------- genes : List[str] @@ -28,6 +36,13 @@ def import_mt( hl.matrixtable.MatrixTable MT of all samples with all variants located within specified genes """ + warnings.warn( + "import_mt() is hard deprecated. hl.import_gvcfs() is no longer available. " + "See WGS 01_QC notebook implementation for alternatives.", + DeprecationWarning, + stacklevel=2, + ) + get_vcfs = partial( lookup_vcfs, mapping=mapping, @@ -65,7 +80,7 @@ def get_position(gene: str, mapping: pd.DataFrame): return chromosome, blocks, start, end -def lookup_regions(gene: str, mapping: pd.DataFrame) -> hl.expr.LocusExpression: +def lookup_regions(gene: str, mapping: pd.DataFrame): chromosome, _, start, end = get_position(gene, mapping) region = [ @@ -92,6 +107,12 @@ def smart_split_multi_mt( mt: hl.matrixtable.MatrixTable, left_aligned=False ) -> hl.matrixtable.MatrixTable: """Split multiple alleles into bi-allelic in a clever way + + .. deprecated:: + This function is soft deprecated and will be removed in a future version. + ML-corrected is already bi-allelic. + + Parameters ---------- mt : hl.matrixtable.MatrixTable @@ -103,6 +124,11 @@ def smart_split_multi_mt( hl.matrixtable.MatrixTable MT with only bi-allelic sites """ + warnings.warn( + "smart_split_multi_mt() is soft deprecated. ML-corrected is already bi-allelic", + DeprecationWarning, + stacklevel=2, + ) mt = mt.key_rows_by("locus", "alleles") @@ -115,3 +141,58 @@ def smart_split_multi_mt( mt = split.union_rows(bi) return mt + + +def make_single_gene_setlist( + mt: hl.matrixtable.MatrixTable, set_list_file: str = "/tmp/GIPR.setlist" +): + """Create a SETLIST file for VEP annotation from a Hail MatrixTable + + .. warning:: + This function only supports single-gene MatrixTables. For multi-gene MTs, + filter to one gene first or call this function separately for each gene. + + Parameters + ---------- + mt : hl.matrixtable.MatrixTable + Hail MatrixTable containing variants from a SINGLE gene + set_list_file : str + Output path for the setlist file + + Returns + ------- + None + Writes a SETLIST file to specified path + + Raises + ------ + ValueError + If the MatrixTable contains variants from multiple genes or chromosomes + """ + # Validate single gene/chromosome + gene_symbols = mt.aggregate_rows( + hl.agg.collect_as_set(mt.vep.transcript_consequences.gene_symbol) + ) + chromosomes = mt.aggregate_rows(hl.agg.collect_as_set(mt.locus.contig)) + + if len(gene_symbols) > 1: + raise ValueError( + f"MatrixTable contains {len(gene_symbols)} genes: {gene_symbols}. " + "make_setlist() only supports single-gene MTs. Filter to one gene first." + ) + + if len(chromosomes) > 1: + raise ValueError( + f"MatrixTable spans {len(chromosomes)} chromosomes: {chromosomes}. " + "make_setlist() requires all variants on the same chromosome." + ) + + # SETLIST file + position = mt.aggregate_rows(hl.agg.min(mt.locus.position)) + names = mt.varid.collect() + names_str = ",".join(names) + + line = f"{list(gene_symbols)[0]}\t{list(chromosomes)[0]}\t{position}\t{names_str}" + + with open(set_list_file, "w") as f: + f.write(line)