diff --git a/build/devenv/devenv.Dockerfile b/build/devenv/devenv.Dockerfile index 53bfd757..9fab629a 100644 --- a/build/devenv/devenv.Dockerfile +++ b/build/devenv/devenv.Dockerfile @@ -95,7 +95,9 @@ RUN apt-get update && \ apt-get install -y --no-install-recommends \ vim \ git \ - htop + htop \ + picard-tools \ + openjdk-8-jre-headless # Install dropbear SSH server and keygen tool COPY --from=dropbear /opt/dropbear/dropbear /usr/bin diff --git a/src/rdds/variant_rank_score/grch38-verification.ipynb b/src/rdds/variant_rank_score/grch38-verification.ipynb new file mode 100644 index 00000000..fa3cffeb --- /dev/null +++ b/src/rdds/variant_rank_score/grch38-verification.ipynb @@ -0,0 +1,601 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "63eb996a-ffb0-47cc-94be-64077f3dc48d", + "metadata": {}, + "source": [ + "# GRCh38 Verification\n", + "\n", + "This notebook is used as a verification tool to make sure MUTACC hg19 solved cases behave identically on\n", + "GRCh38 in raredisease v3.0.\n", + "\n", + "## Input required\n", + "- Subset of MUTACC solved cases from CG, see this [samplesheet](https://github.com/Clinical-Genomics/mivmirvalidation/blob/master/samplesheet.csv)\n", + "- Above grch37 cases run through the mivmir validation pipeline, to annotate variants with mivmir gicam scores (grch37 inferences)\n", + "- Subset of above cases where access to fastq data is available, to run raredisease v3 with mivmir gicam annotation enabled (to generate grch38 inferences)\n", + "\n", + "```\n", + "/rdds/grch38-verification/\n", + " - grch38\n", + " - CASE\n", + " - ...\n", + " - hg19\n", + " - CASE.vcf\n", + " - causatives/\n", + " - CASE.vcf\n", + " - references/\n", + "```\n", + "\n", + "## 1. Convert old MIP cases to hg38" + ] + }, + { + "cell_type": "code", + "execution_count": 164, + "id": "b3f2c0c4-f173-4b43-a7dd-a60b8f90a8a0", + "metadata": {}, + "outputs": [], + "source": [ + "%load_ext autoreload\n", + "%autoreload all -l" + ] + }, + { + "cell_type": "code", + "execution_count": 134, + "id": "073c99ef-8e38-453a-a4c8-acdbf55b13b4", + "metadata": {}, + "outputs": [], + "source": [ + "DATA_DIR='/rdds/grch38-verification'" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "384a5d6f-1f10-4328-b46a-d0421f1905b7", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "NAME=\"Ubuntu\"\n", + "VERSION=\"20.04.6 LTS (Focal Fossa)\"\n", + "ID=ubuntu\n", + "ID_LIKE=debian\n", + "PRETTY_NAME=\"Ubuntu 20.04.6 LTS\"\n", + "VERSION_ID=\"20.04\"\n", + "HOME_URL=\"https://www.ubuntu.com/\"\n", + "SUPPORT_URL=\"https://help.ubuntu.com/\"\n", + "BUG_REPORT_URL=\"https://bugs.launchpad.net/ubuntu/\"\n", + "PRIVACY_POLICY_URL=\"https://www.ubuntu.com/legal/terms-and-policies/privacy-policy\"\n", + "VERSION_CODENAME=focal\n", + "UBUNTU_CODENAME=focal\n" + ] + } + ], + "source": [ + "%%script bash\n", + "cat /etc/os*" + ] + }, + { + "cell_type": "code", + "execution_count": 112, + "id": "ebcfb3fb-a1d9-4c77-ad2d-181289e35443", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "'v1.12.1-35-g6b9b554-dirty'" + ] + }, + "execution_count": 112, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "from rdds.lib.git import git_version; git_version()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "933c6657-a62d-4911-91c1-0c8ec109fb5c", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "%%script bash\n", + "# https://broadinstitute.github.io/picard/\n", + "apt-get update && apt-get install -y picard-tools" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e45501c0-89eb-4820-ab96-cbca04532be8", + "metadata": {}, + "outputs": [], + "source": [ + "# For MUTACC causative variants, remove INFO/END as it's conflicting with coordinates in the new genome ref\n", + "# for csv in `realpath *`; do sed -i 's/END\\=[0123456789]*\\;//g' $csv && sed -i 's/\\#\\#INFO= There's a need for a custom ID.\n", + "\n", + "Furthermore, in the causative VCFs converted to grch38, there's no ID, it's set to missing `.`." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "558f166a-fd98-4ecb-a1c8-22120f3455e4", + "metadata": {}, + "outputs": [], + "source": [ + "# Selected subset of cases for verification\n", + "hg19_grch38_map = {\n", + "'actualgriffon': 'engagedmite',\n", + "'allowedsawfly':'lenientchimp',\n", + "'evidentpony':'lastinggrizzly',\n", + "'factualnewt':'ablecrow',\n", + "'finerram': 'mainantelope',\n", + "'genuinereptile': 'relaxedkiwi',\n", + "'informedmollusk': 'preciousmuskox',\n", + "'nexthare':'humblecow',\n", + "#'patientllama':'learningrooster', Could not retrieve fastq data from CG backend\n", + "'precisespaniel':'amplegrouper'\n", + "}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ec78db09-2ea1-4e35-aa0d-d40bebb35b00", + "metadata": {}, + "outputs": [], + "source": [ + "from glob import glob\n", + "import os\n", + "import pandas as pd\n", + "from tempfile import mkdtemp\n", + "from progressbar import ProgressBar\n", + "import matplotlib.pyplot as plt\n", + "import progressbar.widgets\n", + "import numpy as np\n", + "from cyvcf2 import Writer\n", + "import gc\n", + "\n", + "from rdds.variant_rank_score.inference_exploration.view_vcf_rank_results import view_vcf_rank_results\n", + "from rdds.variant_rank_score.inference_exploration.statfns import plot_performance_vs_threshold\n", + "from rdds.lib.vcf import VCFReader, ParsableVariant" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bdaa7796-a7b3-49c4-b65e-234094d64619", + "metadata": {}, + "outputs": [], + "source": [ + "def vcf_to_pandas(vcf_path: str, is_causative=False) -> pd.DataFrame:\n", + " vcf_reader = VCFReader(vcf_path)\n", + " csq_desc = vcf_reader.csq_description\n", + " n_variants = vcf_reader.number_of_variants\n", + " print(f\"{vcf_path}: {n_variants} variants\")\n", + " chrom = []\n", + " pos = []\n", + " alt = []\n", + " ref = []\n", + " mivmir = []\n", + " gicam = []\n", + " genmod = []\n", + " parse_only_fields = ['POS', 'CHROM', 'RankScoreNormalized', 'MivmirScore', 'GicamScore']\n", + " pbar = ProgressBar(max_value=n_variants)\n", + " for i, variant in enumerate(vcf_reader):\n", + " parsed_variant = ParsableVariant(variant=variant,\n", + " parse_only_fields=parse_only_fields,\n", + " vep_csq_description=vcf_reader.csq_description)\n", + " chrom.append(parsed_variant.CHROM)\n", + " pos.append(int(parsed_variant.POS))\n", + " ref.append(parsed_variant.REF)\n", + " alt.append(parsed_variant.ALT)\n", + " try:\n", + " mivmir.append(parsed_variant.MivmirScore)\n", + " gicam.append(parsed_variant.GicamScore)\n", + " genmod.append(parsed_variant.RankScoreNormalized_value)\n", + " except AttributeError:\n", + " pass\n", + " pbar.update(i)\n", + " pbar.finish()\n", + "\n", + " df = pd.DataFrame({\n", + " 'chrom': chrom,\n", + " 'pos': pos,\n", + " 'ref': ref,\n", + " 'alt': alt,\n", + " 'mivmir': mivmir if len(mivmir) == n_variants else None,\n", + " 'gicam': gicam if len(gicam) == n_variants else None, \n", + " 'genmod': genmod if len(genmod) == n_variants else None,\n", + " 'causative': [1.0] * n_variants if is_causative else [0.0] * n_variants\n", + " })\n", + "\n", + " # Store variant vcf position\n", + " df['vcf_index'] = df.index\n", + " \n", + " # Set index to custom ID\n", + " df.index = df.apply(lambda row: f\"{row.chrom}_{row.pos}-{row.ref}-{row.alt}\", axis=1)\n", + " \n", + " return df\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": 350, + "id": "cfd12a8c-5dfd-40fb-859b-6ea8bafbbf5e", + "metadata": {}, + "outputs": [], + "source": [ + "def _write_variants_to_vcf(indexes: np.ndarray,\n", + " file_name: str,\n", + " output_file_name: str):\n", + " \"\"\"\n", + " Copy variants from input VCF to a new VCF as a subset based on indexes\n", + " \"\"\"\n", + " pbar = ProgressBar(widgets=[progressbar.widgets.BouncingBar()],\n", + " prefix=file_name)\n", + " pbar.start()\n", + " indexes = list(indexes)\n", + " vcf_reader = VCFReader(fname=file_name)\n", + " vcf_writer: Writer = Writer(output_file_name, vcf_reader, mode='w')\n", + " variants = list(vcf_reader) # Load to RAM\n", + " # Write variants in sorted order to output VCF\n", + " for index in indexes:\n", + " vcf_writer.write_record(variants[index])\n", + " vcf_reader.close()\n", + " vcf_writer.close()\n", + " del variants\n", + " gc.collect()\n", + " pbar.finish()" + ] + }, + { + "cell_type": "code", + "execution_count": 368, + "id": "37108e5d-2840-4454-976f-c12e7b105524", + "metadata": {}, + "outputs": [], + "source": [ + "def compute_metrics_per_case(mivmirval_case_name, raredisease_case_name):\n", + " print(f\"hg19 {mivmirval_case_name}:{raredisease_case_name} grch38\")\n", + " mixed_case_name = f\"{mivmirval_case_name}-{raredisease_case_name}\"\n", + " base_work_dir = f'/rdds/tmp/grch38-verification/{mixed_case_name}'\n", + " os.makedirs(base_work_dir, exist_ok=True)\n", + " ref_dir = base_work_dir + '/mivmirvalidation'\n", + " rd_dir = base_work_dir + '/raredisease_v3'\n", + " print(base_work_dir)\n", + " causative_vcf_as_grch38 = os.path.join(DATA_DIR, 'hg19', 'causatives', 'selected-subset', 'as-grch38', f'{mivmirval_case_name}_causative.vcf')\n", + " assert os.path.exists(causative_vcf_as_grch38), causative_vcf_as_grch38\n", + " vcf_mivmir_validation = glob(os.path.join(DATA_DIR, 'hg19', 'as-grch38', mivmirval_case_name + '*.vcf'))[0]\n", + " assert os.path.exists(vcf_mivmir_validation), vcf_mivmir_validation\n", + " vcf_raredisease_grch38 = os.path.join(DATA_DIR, 'grch38', raredisease_case_name, 'rank_and_filter', f'{raredisease_case_name}_snv_ranked_clinical.vcf.gz')\n", + " assert os.path.exists(vcf_raredisease_grch38), vcf_raredisease_grch38\n", + " #view_vcf_rank_results(vcf_file_path=vcf_mivmir_validation, vcf_pathogenic_path=causative_vcf_as_grch38, workdir=ref_dir)\n", + " #view_vcf_rank_results(vcf_file_path=vcf_raredisease_grch38, vcf_pathogenic_path=causative_vcf_as_grch38, workdir=rd_dir)\n", + " df_causative = vcf_to_pandas(causative_vcf_as_grch38, is_causative=True)\n", + " assert len(df_causative) == 1, len(df_causative)\n", + " df_hg19 = vcf_to_pandas(vcf_mivmir_validation)\n", + " df_grch38 = vcf_to_pandas(vcf_raredisease_grch38)\n", + "\n", + " df_causative.to_csv(os.path.join(base_work_dir, 'causative.csv'))\n", + "\n", + " # Store sorted grch38 variants VCF\n", + " store_n_variants = 4000\n", + " for model in ['gicam', 'genmod']:\n", + " _write_variants_to_vcf(indexes=df_grch38.set_index('vcf_index').sort_values(model,\n", + " ascending=False).iloc[0:store_n_variants].index.values,\n", + " file_name=vcf_raredisease_grch38,\n", + " output_file_name=os.path.join(base_work_dir, os.path.basename(vcf_raredisease_grch38).replace('.vcf', f'-grch38-sorted-{model}.vcf')))\n", + "\n", + " # Set causative variant label\n", + " df_hg19.loc[df_causative.index[0], 'causative'] = 1.0\n", + " df_grch38.loc[df_causative.index[0], 'causative'] = 1.0\n", + " assert len(df_hg19[df_hg19.causative == 1]) == 1\n", + " assert len(df_grch38[df_grch38.causative == 1]) == 1\n", + "\n", + " # Performance GRCH38 per model on all of variants in grch38 VCF\n", + " for model in ['gicam', 'mivmir', 'genmod']:\n", + " plot_performance_vs_threshold(predictions=df_grch38[f'{model}'].values,\n", + " labels=df_grch38.causative.values,\n", + " output_path=os.path.join(base_work_dir, f'performance-{model}-all-grch38-variants.png'))\n", + " \n", + " # Merge VCFs based on the intersection of index CHR_POS_REF_ALT in the VCFs (dropping non common variants)\n", + " # This is NOT merging index on different reference genome positions, the intersection refers to variants overlapping MUTACC excerpt and rd v3 fastq inferences\n", + " cols = ['genmod', 'gicam', 'mivmir', 'causative']\n", + " merged_hg19_grch38 = df_hg19[cols].merge(df_grch38[cols], how='inner', left_index=True, right_index=True, suffixes=('_hg19', '_grch38'))\n", + "\n", + " # Compute performance on intersecting variants\n", + " for ref_genome in ['hg19', 'grch38']:\n", + " for model in ['genmod', 'mivmir', 'gicam']:\n", + " plot_performance_vs_threshold(predictions=merged_hg19_grch38[f'{model}_{ref_genome}'].values,\n", + " labels=merged_hg19_grch38[f'causative_{ref_genome}'].values,\n", + " output_path=os.path.join(base_work_dir, f'performance-intersecting-variants-{model}-{ref_genome}.png'))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "df68e974-a23d-49f2-8e03-626a0b6c9884", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "import traceback\n", + "import sys\n", + "\n", + "for mivmirval_case_name, raredisease_case_name in hg19_grch38_map.items():\n", + " try:\n", + " compute_metrics_per_case(mivmirval_case_name, raredisease_case_name)\n", + " except Exception as e:\n", + " print('ERROR', mivmirval_case_name, raredisease_case_name , e)\n", + " traceback.print_exc(file=sys.stdout)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bebd4996-4066-43d9-81d6-2eaf68ffb8f2", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "from glob import glob\n", + "from IPython.display import Image\n", + "genmod_perf_plots = glob('/rdds/tmp/grch38-verification/**/*genmod-all-grch38-variants.png')\n", + "gicam_perf_plots = glob('/rdds/tmp/grch38-verification/**/*gicam-all-grch38-variants.png')\n", + "for genmod_plot, gicam_plot in zip(genmod_perf_plots, gicam_perf_plots):\n", + " print(genmod_plot)\n", + " display(Image(genmod_plot))\n", + " print(gicam_plot)\n", + " display(Image(gicam_plot))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4ed4b6fc-4273-4a7d-aa19-100ec0ab2c19", + "metadata": {}, + "outputs": [], + "source": [ + "from rdds.variant_rank_score.model.model import FEATURES_TEXT, FEATURES_FLOAT\n", + "MODEL_FEATURES = []\n", + "MODEL_FEATURES.extend(FEATURES_TEXT)\n", + "MODEL_FEATURES.extend(FEATURES_FLOAT)\n", + "MODEL_FEATURES.append('GNOMADAF_grpmax') # Replaces GNOMADAF_popmax in grch38\n", + "MODEL_FEATURES.append('MivmirScore') # Store inference to CSV as well\n", + "MODEL_FEATURES.append('GicamScore')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0937c16a-e015-44d5-8385-49fb87d573f9", + "metadata": {}, + "outputs": [], + "source": [ + "# Load causative variants and compare feature inputs\n", + "\n", + "def variant_id(variant) -> str:\n", + " return f\"{variant.CHROM}-{variant.POS}-{variant.REF}-{variant.ALT}\"\n", + "\n", + "for mivmirval_case_name, raredisease_case_name in hg19_grch38_map.items():\n", + " mixed_case_name = f\"{mivmirval_case_name}-{raredisease_case_name}\"\n", + " base_work_dir = f'/rdds/tmp/grch38-verification/{mixed_case_name}'\n", + " print(mixed_case_name)\n", + " df_path = os.path.join(base_work_dir, 'comparison-input-features.csv')\n", + " if os.path.exists(df_path):\n", + " print('Comparison exists, continuing.')\n", + " continue\n", + " causative_file_hg19 = glob(f'/rdds/grch38-verification/hg19/causatives/selected-subset/as-grch38/{mivmirval_case_name}_causative.vcf')[0]\n", + " causative_reader= VCFReader(causative_file_hg19, unpack_if_gzipped=False)\n", + " assert causative_reader.number_of_variants == 1\n", + " causative_variant = list(causative_reader)[0]\n", + " causative_reader.close()\n", + " causative_variant_id = variant_id(causative_variant)\n", + " print('causative:', causative_file_hg19, causative_variant_id)\n", + " # Load mivmirval and raredisease v3 inferences\n", + " print('Reading mivmirval data ...')\n", + " hg19_case_vcf = glob(f'/rdds/grch38-verification/hg19/as-grch38/{mivmirval_case_name}*.vcf')[0]\n", + " print(hg19_case_vcf)\n", + " mivmirval_reader = VCFReader(hg19_case_vcf, unpack_if_gzipped=False)\n", + " mivmirval_causative = None\n", + " for variant in mivmirval_reader:\n", + " if variant_id(variant) == causative_variant_id:\n", + " print('match mivmirval', variant.CHROM, causative_variant_id, variant_id(variant))\n", + " mivmirval_causative = ParsableVariant(variant, vep_csq_description=mivmirval_reader.csq_description)\n", + " break\n", + " assert mivmirval_causative is not None\n", + " mivmirval_reader.close()\n", + " print('Reading raredisease v3 data ...')\n", + " grch38_case_vcf = glob(f'/rdds/grch38-verification/grch38/{raredisease_case_name}/rank_and_filter/{raredisease_case_name}_snv_ranked_clinical.vcf.gz')[0]\n", + " print(grch38_case_vcf)\n", + " raredisease_reader = VCFReader(grch38_case_vcf, unpack_if_gzipped=False)\n", + " raredisease_causative = None\n", + " for variant in raredisease_reader:\n", + " if variant_id(variant) == causative_variant_id:\n", + " print('match RD', variant.CHROM, causative_variant_id, variant_id(variant))\n", + " raredisease_causative = ParsableVariant(variant, vep_csq_description=raredisease_reader.csq_description)\n", + " break\n", + " assert raredisease_causative is not None\n", + " raredisease_reader.close()\n", + "\n", + " df = pd.DataFrame()\n", + " for feature_name in MODEL_FEATURES:\n", + " try:\n", + " mivmireval_feature = getattr(mivmirval_causative, feature_name)\n", + " except AttributeError:\n", + " mivmireval_feature = None\n", + " try:\n", + " raredisease_feature = getattr(raredisease_causative, feature_name)\n", + " except AttributeError:\n", + " raredisease_feature = None\n", + " \n", + " _df = pd.DataFrame(data = {\n", + " feature_name: [mivmireval_feature, raredisease_feature]\n", + " },\n", + " index=['mivmireval-'+variant_id(mivmirval_causative), 'rd-'+variant_id(mivmirval_causative)])\n", + " df = pd.concat((df, _df), axis=1) # concat columnwise \n", + " print('storing', df_path)\n", + " df.to_csv(df_path)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "96095000-ce77-4ad3-b040-65af29cb219b", + "metadata": {}, + "outputs": [], + "source": [ + "df_comparisons = glob('/rdds/tmp/grch38-verification/*/comparison-input-features.csv', recursive=True)\n", + "pd.options.display.max_columns = None\n", + "pd.options.display.max_rows = None\n", + "pd.options.display.max_colwidth = None\n", + "for df_path in df_comparisons:\n", + " print(df_path)\n", + " display(pd.read_csv(df_path))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7fd9755a-09a7-4c0c-b044-690ebfae1024", + "metadata": {}, + "outputs": [], + "source": [ + "# transfer images, preserving catalog names: rsync -v --relative vm:CATALOG_DIR/tmp/grch38-verification/./**/*.png ." + ] + } + ], + "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.8.10" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}