code/Oncoplot.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "fb8d8b61-2443-49cc-9937-67b40cf8d2c0",
   "metadata": {},
   "source": [
    "# Read & parse snpEFF data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed05dca8-b6e6-4a87-9b52-43456d949777",
   "metadata": {},
   "outputs": [],
   "source": [
    "_SAMPLE_TABLE = pd.read_excel(\"<sampletable\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b5cbc22-a71e-4e53-b1c9-a53be0fabf2e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import pathlib\n",
    "import tqdm.auto as tqdm\n",
    "\n",
    "tqdm.tqdm.pandas()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b2b4ee8-a945-43e4-a2fd-fe4a8f3427a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "files = pd.Series(\n",
    "    pathlib.Path(\"variantcalling/annotation/\").rglob(\"*snpEff.ann.vcf.gz\")\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46ec0f2e-7ccd-4019-8f59-cc9f09f77f21",
   "metadata": {},
   "outputs": [],
   "source": [
    "files.index = files.map(lambda x: x.parent.name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f1bfe9d-4132-45d1-b250-969b67677744",
   "metadata": {},
   "outputs": [],
   "source": [
    "genes = [\n",
    "    \"LZTR1\",\n",
    "    \"SMARCA4\",\n",
    "    \"PDGFRA\",\n",
    "    \"LDB1\",\n",
    "    \"EIF1AX\",\n",
    "    \"MSH6\",\n",
    "    \"HDAC2\",\n",
    "    \"PTPMT1\",\n",
    "    \"TRAF7\",\n",
    "    \"ALK\",\n",
    "    \"PTPN11\",\n",
    "    \"NF1\",\n",
    "    \"SMARCE1\",\n",
    "    \"HRAS\",\n",
    "    \"PTCH1\",\n",
    "    \"KIT\",\n",
    "    \"ATRX\",\n",
    "    \"AIFM3\",\n",
    "    \"PRKAR1A\",\n",
    "    \"FLT3\",\n",
    "    \"NDUFS3\",\n",
    "    \"SPRED1\",\n",
    "    \"ACVR1\",\n",
    "    \"BRPF1\",\n",
    "    \"D2HGDH\",\n",
    "    \"MYB\",\n",
    "    \"PALB2\",\n",
    "    \"CHD7\",\n",
    "    \"MET\",\n",
    "    \"DGCR8\",\n",
    "    \"KLF4\",\n",
    "    \"YAP1\",\n",
    "    \"CHEK2\",\n",
    "    \"PHOX2B\",\n",
    "    \"MLL2\",\n",
    "    \"PPL\",\n",
    "    \"RAF1\",\n",
    "    \"SMARCD1\",\n",
    "    \"MRE11A\",\n",
    "    \"ASXL1\",\n",
    "    \"NF2\",\n",
    "    \"NBN\",\n",
    "    \"CIC\",\n",
    "    \"ATXN1L\",\n",
    "    \"SDCCAG8\",\n",
    "    \"SUZ12\",\n",
    "    \"HTRA1\",\n",
    "    \"CDKN2C\",\n",
    "    \"MRE11\",\n",
    "    \"DROSHA\",\n",
    "    \"BRCA2\",\n",
    "    \"CDKN2A\",\n",
    "    \"PIK3C2B\",\n",
    "    \"CDK6\",\n",
    "    \"STAT6\",\n",
    "    \"NRAS\",\n",
    "    \"ZMYM3\",\n",
    "    \"IST1\",\n",
    "    \"EPAS1\",\n",
    "    \"PIK3CA\",\n",
    "    \"KDM6A\",\n",
    "    \"PAX6\",\n",
    "    \"GNAQ\",\n",
    "    \"HIST1H3B\",\n",
    "    \"IKBKAP\",\n",
    "    \"MYC\",\n",
    "    \"BRCA1\",\n",
    "    \"PCDH8\",\n",
    "    \"TP53\",\n",
    "    \"DAXX\",\n",
    "    \"PIK3C2G\",\n",
    "    \"H3F3A\",\n",
    "    \"HIST1H3C\",\n",
    "    \"COL2A1\",\n",
    "    \"PIK3R2\",\n",
    "    \"AIP\",\n",
    "    \"SH3PXD2A\",\n",
    "    \"ARID1B\",\n",
    "    \"PLCH1\",\n",
    "    \"MSH2\",\n",
    "    \"MN1\",\n",
    "    \"POLE\",\n",
    "    \"GNA11\",\n",
    "    \"PIK3CG\",\n",
    "    \"SMARCD2\",\n",
    "    \"AKT2\",\n",
    "    \"KMT2B\",\n",
    "    \"BCOR\",\n",
    "    \"GNAS\",\n",
    "    \"KMT2C\",\n",
    "    \"IDH2\",\n",
    "    \"CDH1\",\n",
    "    \"PLXNB1\",\n",
    "    \"STAG2\",\n",
    "    \"FGFR2\",\n",
    "    \"FGFR1\",\n",
    "    \"FOXO3\",\n",
    "    \"KBTBD4\",\n",
    "    \"MYCN\",\n",
    "    \"MGMT\",\n",
    "    \"SMAD4\",\n",
    "    \"PIK3R1\",\n",
    "    \"ZIC1\",\n",
    "    \"H3F3B\",\n",
    "    \"H2AX\",\n",
    "    \"PTCH2\",\n",
    "    \"APC\",\n",
    "    \"FGFR3\",\n",
    "    \"MDM2\",\n",
    "    \"MEN1\",\n",
    "    \"ATM\",\n",
    "    \"TSPAN31\",\n",
    "    \"TERT\",\n",
    "    \"FGFR4\",\n",
    "    \"ARID2\",\n",
    "    \"ARID1A\",\n",
    "    \"SMARCB1\",\n",
    "    \"NAB2\",\n",
    "    \"KMT2D\",\n",
    "    \"H3P6\",\n",
    "    \"LRRFIP2\",\n",
    "    \"MYBL1\",\n",
    "    \"RET\",\n",
    "    \"PMS2\",\n",
    "    \"JAK3\",\n",
    "    \"CCND2\",\n",
    "    \"CSF1R\",\n",
    "    \"H3C3\",\n",
    "    \"DICER1\",\n",
    "    \"MLH1\",\n",
    "    \"KLK1\",\n",
    "    \"IDH1\",\n",
    "    \"ZIC4\",\n",
    "    \"CTNNB1\",\n",
    "    \"RAD50\",\n",
    "    \"TBR1\",\n",
    "    \"BRAF\",\n",
    "    \"KCNK12\",\n",
    "    \"SMARCD3\",\n",
    "    \"ERBB2\",\n",
    "    \"ROS1\",\n",
    "    \"FOXR2\",\n",
    "    \"IDO2\",\n",
    "    \"EZH2\",\n",
    "    \"JAK2\",\n",
    "    \"MDM4\",\n",
    "    \"EZHIP\",\n",
    "    \"MPL\",\n",
    "    \"MYL1\",\n",
    "    \"NOTCH1\",\n",
    "    \"ATXN1\",\n",
    "    \"C11ORF95\",\n",
    "    \"GSE1\",\n",
    "    \"CDK4\",\n",
    "    \"GABRA6\",\n",
    "    \"BRPF3\",\n",
    "    \"BAP1\",\n",
    "    \"EED\",\n",
    "    \"EGFR\",\n",
    "    \"HNF1A\",\n",
    "    \"DDX3X\",\n",
    "    \"RBPJ\",\n",
    "    \"TH2LCRR\",\n",
    "    \"CREBBP\",\n",
    "    \"MYCNOS\",\n",
    "    \"ABL1\",\n",
    "    \"CCND1\",\n",
    "    \"SMARCA2\",\n",
    "    \"SUFU\",\n",
    "    \"TSC1\",\n",
    "    \"TSC2\",\n",
    "    \"SMIM4\",\n",
    "    \"SETD2\",\n",
    "    \"H2AFX\",\n",
    "    \"PBRM1\",\n",
    "    \"SF3B1\",\n",
    "    \"DDR1\",\n",
    "    \"FBXW7\",\n",
    "    \"AKT3\",\n",
    "    \"TCF4\",\n",
    "    \"PRKCA\",\n",
    "    \"KDR\",\n",
    "    \"H3C2\",\n",
    "    \"SMO\",\n",
    "    \"KRAS\",\n",
    "    \"FUBP1\",\n",
    "    \"NTRK2\",\n",
    "    \"RB1\",\n",
    "    \"AKT1\",\n",
    "    \"NDRG2\",\n",
    "    \"PTEN\",\n",
    "    \"ESR1\",\n",
    "    \"ATR\",\n",
    "    \"GPR161\",\n",
    "    \"SRP19\",\n",
    "    \"VHL\",\n",
    "    \"PPM1D\",\n",
    "    \"RELA\",\n",
    "    \"CTDNEP1\",\n",
    "    \"MTOR\",\n",
    "    \"POLD1\",\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e83c529-189c-45f6-abf1-476b97e5a7b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_and_filter(file):\n",
    "    tmp = pd.read_csv(file, comment=\"#\", sep=\"\\t\", header=None)\n",
    "    tmp = tmp[tmp[6] == \"PASS\"].drop(columns=[2, 6, 8])\n",
    "    tmp2 = pd.DataFrame()\n",
    "    if not tmp.empty:\n",
    "        tmp2 = tmp[\n",
    "            (tmp[7].str.split(\"|\", expand=True)[3] == \"TERT\")\n",
    "            & (tmp[7].str.split(\"|\", expand=True)[1] == \"upstream_gene_variant\")\n",
    "        ]\n",
    "        tmp = tmp[tmp[7].str.split(\"|\", expand=True)[2].isin([\"HIGH\", \"MODERATE\"])]\n",
    "    if not tmp.empty:\n",
    "        genes_idx = (\n",
    "            tmp[7]\n",
    "            .str.split(\"|\", expand=True)[3]\n",
    "            .str.split(\"&\", expand=True)\n",
    "            .apply(lambda x: x.isin(genes))\n",
    "            .any(axis=1)\n",
    "        )\n",
    "        tmp = tmp[genes_idx]\n",
    "    tmp = pd.concat((tmp2, tmp)).drop_duplicates()\n",
    "    return tmp.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58e798d6-a5ff-45fc-9404-25e869a3c28b",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = files.progress_apply(read_and_filter)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a26f369e-2572-4a14-af10-25dae3846343",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.concat(data.to_dict())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa54ff68-c2e5-42ed-92cd-4cb71b9b6f5e",
   "metadata": {},
   "outputs": [],
   "source": [
    "data[\"Gene\"] = data[7].str.split(\"|\", expand=True)[3]\n",
    "data[\"Impact\"] = data[7].str.split(\"|\", expand=True)[2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "944bf473-b1a9-4aac-bbb4-6fbde7678a84",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = data.reset_index(level=1, drop=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12b10f75-54af-4a6a-91b7-5a4cf8f30955",
   "metadata": {},
   "outputs": [],
   "source": [
    "data[\"variant\"] = data[7].str.split(\"|\", expand=True)[9]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e61c75d-9450-49fa-9c64-f88bb5f76fee",
   "metadata": {},
   "outputs": [],
   "source": [
    "data.to_csv(\"mutation_status_snpEff.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2491754f-98be-4417-a589-d1af2ba2baec",
   "metadata": {},
   "source": [
    "# Read & parse VEP data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fdc68200-d46d-4be2-bcf6-e5aaeb9ade67",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import pathlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8c8806f-224d-49c3-918e-157fbd3e0b41",
   "metadata": {},
   "outputs": [],
   "source": [
    "files = pd.Series(pathlib.Path(\"variantcalling/filtered/\").glob(\"*VEP*\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87813f85-236b-4744-b5dd-eb55263b628a",
   "metadata": {},
   "outputs": [],
   "source": [
    "files.index = files.map(lambda x: x.name).str.split(\".\", n=1, expand=True)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1eff5034-de39-4df4-936d-eef208c481e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "cols = \"Allele|Consequence|IMPACT|SYMBOL|Gene|Feature_type|Feature|BIOTYPE|EXON|INTRON|HGVSc|HGVSp|cDNA_position|CDS_position|Protein_position|Amino_acids|Codons|Existing_variation|DISTANCE|STRAND|FLAGS|VARIANT_CLASS|SYMBOL_SOURCE|HGNC_ID|CANONICAL|MANE_SELECT|MANE_PLUS_CLINICAL|TSL|APPRIS|CCDS|ENSP|SWISSPROT|TREMBL|UNIPARC|UNIPROT_ISOFORM|GENE_PHENO|SIFT|PolyPhen|DOMAINS|miRNA|AF|AFR_AF|AMR_AF|EAS_AF|EUR_AF|SAS_AF|gnomADe_AF|gnomADe_AFR_AF|gnomADe_AMR_AF|gnomADe_ASJ_AF|gnomADe_EAS_AF|gnomADe_FIN_AF|gnomADe_NFE_AF|gnomADe_OTH_AF|gnomADe_SAS_AF|gnomADg_AF|gnomADg_AFR_AF|gnomADg_AMI_AF|gnomADg_AMR_AF|gnomADg_ASJ_AF|gnomADg_EAS_AF|gnomADg_FIN_AF|gnomADg_MID_AF|gnomADg_NFE_AF|gnomADg_OTH_AF|gnomADg_SAS_AF|MAX_AF|MAX_AF_POPS|FREQS|CLIN_SIG|SOMATIC|PHENO|PUBMED|MOTIF_NAME|MOTIF_POS|HIGH_INF_POS|MOTIF_SCORE_CHANGE|TRANSCRIPTION_FACTORS\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4bfb8d2d-1a02-4c10-b63d-58a18a37a878",
   "metadata": {},
   "outputs": [],
   "source": [
    "cols = cols.split(\"|\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b80d3dd-a66b-4c66-9763-97d803198a7e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def read(file):\n",
    "    try:\n",
    "        tmp = pd.read_csv(\n",
    "            file,\n",
    "            comment=\"#\",\n",
    "            sep=\"\\t\",\n",
    "            header=None,\n",
    "            index_col=[0, 1, 2, 3, 4],\n",
    "            usecols=[0, 1, 2, 3, 4, 7],\n",
    "        )[7]\n",
    "    except pd.errors.EmptyDataError:\n",
    "        return pd.DataFrame()\n",
    "    tmp = tmp.str.extract(\"CSQ=([^;]+)\")\n",
    "    tmp = tmp[0].str.split(\"|\", expand=True, n=len(cols))\n",
    "    if len(tmp.columns) == 1:\n",
    "        return pd.DataFrame()\n",
    "    tmp.columns = (cols + [\"other\"]) if len(tmp.columns) > len(cols) else cols\n",
    "    return tmp[tmp.IMPACT.isin([\"HIGH\", \"MODERATE\"])].copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46878d37-48e6-489a-a68a-9f895394b1fe",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = files.map(read)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70d1f0db-3e1f-4fcc-865c-9d7d1ba7b759",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.concat(data.to_dict())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "43b3f955-dae6-4855-9370-708e3494e85e",
   "metadata": {},
   "outputs": [],
   "source": [
    "data.to_csv(\"vep_annotation.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bfddfd9b-9b66-44c0-8be1-4ade5a641453",
   "metadata": {},
   "source": [
    "# Plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd2e5065-75d8-42bf-8bbe-6ebf55cddfd1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae39f297-3c03-4790-b73c-97fff0217a55",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = _SAMPLE_TABLE.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b10e1d75-8626-41cb-8f54-cc8af8ffca57",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = data.set_index(\"Supplier_Nummer\")[\n",
    "    [\n",
    "        \"Age\",\n",
    "        \"Neural Subgroup\",\n",
    "        \"Sex\",\n",
    "        \"DKFZ methylation subclass\",\n",
    "        \"MGMT\",\n",
    "        \"Location\",\n",
    "        \"Idat\",\n",
    "    ]\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6d25373-7bd8-4765-943b-42e18b07d32c",
   "metadata": {},
   "outputs": [],
   "source": [
    "palette_group = {\n",
    "    \"Sex\": {\n",
    "        \"male\": sns.color_palette(\"tab20\")[0],\n",
    "        \"female\": sns.color_palette(\"tab20\")[2],\n",
    "    },\n",
    "    \"DKFZ methylation subclass\": {\n",
    "        \"RTK 1\": sns.color_palette(\"tab20\")[8],\n",
    "        \"RTK 2\": sns.color_palette(\"tab20\")[9],\n",
    "        \"MES\": sns.color_palette(\"tab20\")[7],\n",
    "    },\n",
    "    \"MGMT\": {\n",
    "        \"methylated\": sns.color_palette(\"Greys\")[5],\n",
    "        \"non-methylated\": sns.color_palette(\"Greys\")[0],\n",
    "    },\n",
    "    \"Location\": {\n",
    "        \"temporal\": sns.color_palette(\"tab20\")[10],\n",
    "        \"parietal\": sns.color_palette(\"tab20\")[12],\n",
    "        \"occipital\": sns.color_palette(\"tab20\")[14],\n",
    "        \"frontal\": sns.color_palette(\"tab20\")[16],\n",
    "        \"brain, nos\": sns.color_palette(\"tab20\")[18],\n",
    "        \"cerebrum\": sns.color_palette(\"tab20\")[19],\n",
    "    },\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1239f548-834a-4e74-91c5-2cdf0f21d4e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "palette = {k: v for x in palette_group.values() for k, v in x.items()}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "297c460d-f321-46d9-b068-48ef243442eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set_context(\"talk\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e65fe570-7908-466e-8375-b2dea7c969e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "anno = data.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8288540-37e9-4135-aced-baf46aca49a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_csv(\"mutation_status_snpEff.csv\", index_col=0)\n",
    "\n",
    "data = data[data[\"7\"].str.extract(\"MQ=(\\d+)\")[0].astype(int) > 30]\n",
    "\n",
    "data[\"type\"] = data[\"7\"].str.split(\"|\", expand=True)[1]\n",
    "\n",
    "data[\"type\"].unique()\n",
    "\n",
    "data = data[data[\"5\"] > 20]\n",
    "data = data[pd.to_numeric(data[\"9\"].str.split(\":\", expand=True)[3]) > 15]\n",
    "\n",
    "data[\"consensus_type\"] = pd.NA\n",
    "data = data[\n",
    "    ~data.index.isin(\n",
    "        data.query(\"Gene == 'TERT' & type != 'upstream_gene_variant'\").index\n",
    "    )\n",
    "]\n",
    "data.loc[\n",
    "    data.consensus_type.isna() & (data.variant.str.contains(\"c.-\")), \"consensus_type\"\n",
    "] = \"Promoter\"\n",
    "data.loc[\n",
    "    data.consensus_type.isna() & data.type.str.contains(\"missense_variant\"),\n",
    "    \"consensus_type\",\n",
    "] = \"Missense\"\n",
    "data.loc[\n",
    "    data.consensus_type.isna() & data.type.str.contains(\"stop_gained\"), \"consensus_type\"\n",
    "] = \"Truncation\"\n",
    "data.loc[\n",
    "    data.consensus_type.isna()\n",
    "    & data.variant.str.contains(\"dup\")\n",
    "    & ~data.type.str.contains(\"intron_variant\"),\n",
    "    \"consensus_type\",\n",
    "] = \"InDel\"\n",
    "data.loc[\n",
    "    data.consensus_type.isna()\n",
    "    & data.variant.str.contains(\"del\")\n",
    "    & ~data.type.str.contains(\"intron_variant\"),\n",
    "    \"consensus_type\",\n",
    "] = \"InDel\"\n",
    "data.loc[\n",
    "    data.consensus_type.isna()\n",
    "    & data.variant.str.contains(\"ins\")\n",
    "    & ~data.type.str.contains(\"intron_variant\"),\n",
    "    \"consensus_type\",\n",
    "] = \"InDel\"\n",
    "data.loc[\n",
    "    data.consensus_type.isna() & data.type.str.contains(\"stop_lost\"), \"consensus_type\"\n",
    "] = \"Missense\"\n",
    "data.loc[data.consensus_type.isna(), \"consensus_type\"] = \"Missense\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1d3e50f-c0a7-43c9-987c-5fc326d55862",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = data[~data.type.str.contains(\"intron_variant\")]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0122f554-c693-4071-bbcb-fd4126bcc35d",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = data.set_index([\"0\", \"1\", \"3\", \"4\"], append=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "777bd721-f410-4b4e-9f6f-03a5e0b8a3da",
   "metadata": {},
   "outputs": [],
   "source": [
    "annot = pd.read_csv(\"variantcalling/vep_annotation.csv\", index_col=[0, 1, 2, 4, 5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77489d10-d38d-433a-a14f-1d5c544eba95",
   "metadata": {},
   "outputs": [],
   "source": [
    "data[\"1000G_AF\"] = annot[\"AF\"].reindex(data.index).fillna(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e082a9fe-babf-42e0-b31c-1297a518a646",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = data[~data.variant.isin([\"c.5885C>T\", \"c.2447dupA\"])]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "566eea27-bb1a-49ca-8011-9f4c5463d1d1",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = data.reset_index(level=[1, 2, 3, 4])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "139db2c6-2442-44ce-bd0e-94b921517641",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = data.query(\"`1000G_AF` < 0.001\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8091e60-e016-48c2-b354-cc16f8cfce7a",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = data[~(data.variant == \"c.5885C>T\")]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8431cc0-0ac0-4ebe-99c9-aebefec7c51d",
   "metadata": {},
   "outputs": [],
   "source": [
    "data[\"Supplier_Nummer\"] = (\n",
    "    pd.read_csv(\"ID_mapping.csv\")\n",
    "    .set_index(\"Supplier_Nummer\")[\"Supplier_Nummer.1\"]\n",
    "    .squeeze()\n",
    "    .reindex(data.index)\n",
    "    .values\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14ad8247-786b-4f34-82f4-dc59c0a398f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = data[data.Supplier_Nummer.isin(anno.index)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ddc7eb94-4b31-417b-a72e-3b47859e6a79",
   "metadata": {},
   "outputs": [],
   "source": [
    "mapped = data.pivot_table(\n",
    "    index=\"Supplier_Nummer\", columns=\"Gene\", values=\"consensus_type\", aggfunc=list\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb714c40-9f10-4e8d-9a2d-1dd8e5a623a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "def assign_new(x):\n",
    "    try:\n",
    "        x = list(set(x))\n",
    "        l = len(x)\n",
    "    except:\n",
    "        return pd.NA\n",
    "    if l == 1:\n",
    "        return x[0]\n",
    "    return \"Multiple\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b901b131-dad5-4f76-a2cf-a9601fe2036d",
   "metadata": {},
   "outputs": [],
   "source": [
    "mutations = mapped.applymap(assign_new)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9474f608-437d-4fd0-a4de-01e2178d42bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "mutations = mutations.stack().reset_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d3b3146-4259-43c0-9111-9d1a89ce607a",
   "metadata": {},
   "outputs": [],
   "source": [
    "mutations = mutations.rename(columns={0: \"type\"})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d77e9192-c1c6-4662-a407-130699a02657",
   "metadata": {},
   "outputs": [],
   "source": [
    "genes = [\n",
    "    \"PDGFRA\",\n",
    "    \"PIK3CA\",\n",
    "    \"MDM2\",\n",
    "    \"PTPN11\",\n",
    "    \"NF1\",\n",
    "    \"ATM\",\n",
    "    \"TERT\",\n",
    "    \"TP53\",\n",
    "    \"ATRX\",\n",
    "    \"KIT\",\n",
    "    \"CDK4\",\n",
    "    \"EGFR\",\n",
    "    \"TSC2\",\n",
    "    \"MET\",\n",
    "    \"SETD2\",\n",
    "    \"CCND2\",\n",
    "    \"KDR\",\n",
    "    \"IDH1\",\n",
    "    \"BCOR\",\n",
    "    \"KMT2C\",\n",
    "    \"RB1\",\n",
    "    \"PTEN\",\n",
    "    \"BRAF\",\n",
    "    \"CDKN2A\",\n",
    "    \"FGFR1\",\n",
    "    \"MYCN\",\n",
    "    \"PIK3R1\",\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc6f33d0-6490-447c-aba3-4ae1ee7c592a",
   "metadata": {},
   "outputs": [],
   "source": [
    "mut_select = mutations[mutations.Gene.isin(genes)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94a4978a-0de9-4a76-a137-0a867e525356",
   "metadata": {},
   "outputs": [],
   "source": [
    "mut_select.Gene = pd.Categorical(mut_select.Gene, categories=genes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d6e5f44-5cca-498d-b6d8-0ff4ec39b741",
   "metadata": {},
   "outputs": [],
   "source": [
    "mutfraq = mutations.Supplier_Nummer.value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31ebcbf8-a5ff-4289-8b58-7b82ce239605",
   "metadata": {},
   "outputs": [],
   "source": [
    "mut_palette = sns.color_palette(n_colors=mutations.type.nunique())\n",
    "mut_palette = dict(zip(mutations.type.unique(), mut_palette))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "804e3984-2f97-4435-9125-aa725ceb4457",
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set_context(\"paper\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bcec116b-43b2-49e5-9a6e-e16dacf9760e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_plot(data, mutations, mutfreq, pos=((0, -9), (0, -1.5))):\n",
    "    fig, ax = plt.subplots(\n",
    "        4,\n",
    "        2,\n",
    "        sharex=\"col\",\n",
    "        sharey=\"row\",\n",
    "        gridspec_kw=dict(\n",
    "            width_ratios=[1, 0.1],\n",
    "            height_ratios=[0.5, 2, 0.5, 0.5],  # 0.5],\n",
    "            hspace=0.03,\n",
    "            wspace=0.01,\n",
    "        ),\n",
    "        figsize=(11, 9),\n",
    "    )\n",
    "    axes_dict = {\n",
    "        \"anno\": ax[0, 0],\n",
    "        \"mutation\": ax[1, 0],\n",
    "        \"freq\": ax[1, 1],\n",
    "        \"age\": ax[2, 0],\n",
    "        #  \"CNV\": ax[3,0],\n",
    "        \"patfreq\": ax[3, 0],\n",
    "    }\n",
    "    ax[0, 1].axis(\"off\")\n",
    "    ax[2, 1].axis(\"off\")\n",
    "    ax[3, 1].axis(\"off\")\n",
    "    # ax[4,1].axis(\"off\")\n",
    "    sns.scatterplot(x=data.index, y=data.Age, ax=axes_dict[\"age\"], color=\"lightgray\")\n",
    "    axes_dict[\"age\"].set_xticklabels(\"\")\n",
    "    axes_dict[\"age\"].set_xlabel(\"\")\n",
    "    axes_dict[\"age\"].tick_params(bottom=False)\n",
    "    axes_dict[\"age\"].set_ylabel(\"Age\\n(years)\")\n",
    "    axes_dict[\"age\"].set_ylim(15, 86)\n",
    "    anno = (\n",
    "        data[[\"Sex\", \"DKFZ methylation subclass\", \"MGMT\", \"Location\"]]\n",
    "        .reset_index()\n",
    "        .melt(id_vars=[\"Supplier_Nummer\"])\n",
    "    )\n",
    "    sns.scatterplot(\n",
    "        x=\"Supplier_Nummer\",\n",
    "        y=\"variable\",\n",
    "        ax=axes_dict[\"anno\"],\n",
    "        data=anno,\n",
    "        hue=\"value\",\n",
    "        marker=\"s\",\n",
    "        palette=palette,\n",
    "        legend=True,\n",
    "        s=50,\n",
    "    )\n",
    "    sns.move_legend(axes_dict[\"anno\"], loc=pos[0], ncol=5, title=\"Clinical annotation\")\n",
    "    axes_dict[\"anno\"].set_xticklabels(\"\")\n",
    "    axes_dict[\"anno\"].set_xlabel(\"\")\n",
    "    axes_dict[\"anno\"].set_ylabel(\"\")\n",
    "    axes_dict[\"anno\"].set_ylim(-1, 4)\n",
    "    axes_dict[\"anno\"].tick_params(bottom=False)\n",
    "    # axes_dict[\"anno\"].legend(bbox_to_anchor=axes_dict[\"patfreq\"].bbox)\n",
    "    sns.scatterplot(\n",
    "        x=\"Supplier_Nummer\",\n",
    "        y=\"Gene\",\n",
    "        ax=axes_dict[\"mutation\"],\n",
    "        data=mutations,\n",
    "        hue=\"type\",\n",
    "        marker=\"s\",\n",
    "        legend=True,\n",
    "        palette=mut_palette,\n",
    "        s=50,\n",
    "    )\n",
    "    axes_dict[\"mutation\"].tick_params(bottom=False)\n",
    "    axes_dict[\"mutation\"].set_ylabel(\"\")\n",
    "    sns.move_legend(axes_dict[\"mutation\"], loc=pos[1], ncol=5, title=\"Alteration\")\n",
    "\n",
    "    counts = (\n",
    "        mutations.Gene.value_counts(sort=False) / data.index.size * 100\n",
    "    ).reset_index()\n",
    "    counts.columns = [\"Gene\", \"count\"]\n",
    "    axes_dict[\"freq\"].barh(counts[\"Gene\"], counts[\"count\"], color=\"black\")\n",
    "    axes_dict[\"freq\"].tick_params(left=False)\n",
    "    axes_dict[\"freq\"].set_xticks([0, 10, 20, 30, 40, 50])\n",
    "    axes_dict[\"freq\"].set_xticklabels(\n",
    "        [\"0\\%\", \"10\\%\", \"20\\%\", \"30\\%\", \"40\\%\", \"50\\%\"], rotation=-90\n",
    "    )\n",
    "    axes_dict[\"freq\"].tick_params(labelbottom=True)\n",
    "    mutfreq.name = \"count\"\n",
    "    sns.barplot(\n",
    "        x=\"Supplier_Nummer\",\n",
    "        y=\"count\",\n",
    "        data=mutfreq.reset_index(),\n",
    "        color=\"black\",\n",
    "        ax=axes_dict[\"patfreq\"],\n",
    "    )\n",
    "    axes_dict[\"patfreq\"].set_xticklabels(\n",
    "        axes_dict[\"patfreq\"].get_xticklabels(), rotation=90\n",
    "    )\n",
    "    axes_dict[\"patfreq\"].set_xlabel(\"\")\n",
    "    axes_dict[\"patfreq\"].set_ylabel(\"N. of\\nMutations\")\n",
    "    axes_dict[\"patfreq\"].set_ylim(-0, 73)\n",
    "    return fig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8583ed2b-05c8-469a-8311-f1b6007167b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "subset = anno.query(\"`Neural Subgroup` == 'low'\")\n",
    "mut_subset = mut_select[mut_select.Supplier_Nummer.isin(subset.index)]\n",
    "mutfraq_subset = mutfraq.reindex(subset.index)\n",
    "fig = compute_plot(subset, mut_subset, mutfraq_subset, pos=((0, -10.2), (0, -1.8)))\n",
    "fig.suptitle(\"Neural low\")\n",
    "fig.savefig(\"oncopanel_neural_low.svg\", dpi=200, bbox_inches=\"tight\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a606dd97-2413-4a65-8288-0be137bbbca1",
   "metadata": {},
   "outputs": [],
   "source": [
    "subset = anno.query(\"`Neural Subgroup` == 'high'\")\n",
    "mut_subset = mut_select[mut_select.Supplier_Nummer.isin(subset.index)]\n",
    "mutfraq_subset = mutfraq.reindex(subset.index)\n",
    "fig = compute_plot(subset, mut_subset, mutfraq_subset, pos=((0, -10.2), (0, -1.8)))\n",
    "fig.suptitle(\"Neural high\")\n",
    "fig.savefig(\"oncopanel_neural_high.svg\", dpi=200, bbox_inches=\"tight\")"
   ]
  }
 ],
 "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.9.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}