scripts/extract-chr21-genes-variants.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import gzip"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Extract variants in chr21 into TSV\n",
    "Read VCF and extract coordinates, SV type and variant ID\n",
    "\n",
    "## GCRh38"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "dx download dbVar_byAssembly_hg38/GRCh38.variant_call.all.vcf.gz*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "outf = open('all_variants_chr21.tsv', 'w')\n",
    "# header\n",
    "outf.write('chr\\tstart\\tend\\ttype\\tvariant_id\\n')\n",
    "# read VCF with all variants\n",
    "svtypes = {}\n",
    "for line in gzip.open('GRCh38.variant_call.all.vcf.gz', 'rb'):\n",
    "    line = line.decode('ascii')\n",
    "    # skip headers\n",
    "    if line[0] == '#':\n",
    "        continue\n",
    "    line = line.rstrip().split('\\t')\n",
    "    # skip if not chr 21\n",
    "    if line[0] != '21':\n",
    "        continue\n",
    "    # parse INFO field\n",
    "    infos = {}\n",
    "    for info in line[7].split(';'):\n",
    "        info_pair = info.split('=')\n",
    "        if len(info_pair) == 1:\n",
    "            infos[info] = True\n",
    "        else:\n",
    "            infos[info_pair[0]] = info_pair[1]\n",
    "    # prepare info for TSV\n",
    "    out_line = [line[0]]\n",
    "    # start pos\n",
    "    pos_s = line[1]\n",
    "    out_line.append(pos_s)\n",
    "    # end pos\n",
    "    pos_e = pos_s\n",
    "    if 'END' in infos:\n",
    "        pos_e = infos['END']\n",
    "    out_line.append(pos_e)\n",
    "    # SV type\n",
    "    svtype = infos['SVTYPE']\n",
    "    out_line.append(svtype)\n",
    "    svtypes[svtype] = True\n",
    "    # skip if BND type\n",
    "    if svtype == 'BND':\n",
    "        continue\n",
    "    if int(pos_e) < int(pos_s):\n",
    "        print(svtype)\n",
    "    # variant id\n",
    "    out_line.append(line[2])\n",
    "    # write line in tsv\n",
    "    outf.write('\\t'.join(out_line) + '\\n')\n",
    "outf.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dict_keys(['INS', 'DEL', 'DUP', 'CNV', 'INV', 'BND'])\n"
     ]
    }
   ],
   "source": [
    "print(svtypes.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ID                  file-Fy33VG00Z5b26bpZJ3Zk8f8k\n",
      "Class               file\n",
      "Project             project-Fy1b7V80Z5b4jXb224P1fY4b\n",
      "Folder              /chr-subset-genes-variants\n",
      "Name                all_variants_chr21.tsv\n",
      "State               closing\n",
      "Visibility          visible\n",
      "Types               -\n",
      "Properties          -\n",
      "Tags                -\n",
      "Outgoing links      -\n",
      "Created             Tue Oct 13 23:30:56 2020\n",
      "Created by          jmonlong\n",
      " via the job        job-Fy332xj0Z5bKZz0F8v9gKVfJ\n",
      "Last modified       Tue Oct 13 23:30:57 2020\n",
      "Media type          \n",
      "archivalState       \"live\"\n",
      "cloudAccount        \"cloudaccount-dnanexus\"\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "dx upload all_variants_chr21.tsv --path chr-subset-genes-variants/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GRCh37\n",
    "\n",
    "To match gnomAD-SV frequencies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Error: path \"/opt/notebooks/GRCh37.variant_call.all.vcf.gz\" already exists but\n",
      "-f/--overwrite was not set\n"
     ]
    },
    {
     "ename": "CalledProcessError",
     "evalue": "Command 'b'dx download dbVar_byAssembly_GRCh37/GRCh37.variant_call.all.vcf.gz\\n'' returned non-zero exit status 1.",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mCalledProcessError\u001b[0m                        Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-10-17221335374b>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mget_ipython\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrun_cell_magic\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'bash'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m''\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'dx download dbVar_byAssembly_GRCh37/GRCh37.variant_call.all.vcf.gz\\n'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m/opt/conda/lib/python3.6/site-packages/IPython/core/interactiveshell.py\u001b[0m in \u001b[0;36mrun_cell_magic\u001b[0;34m(self, magic_name, line, cell)\u001b[0m\n\u001b[1;32m   2369\u001b[0m             \u001b[0;32mwith\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbuiltin_trap\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2370\u001b[0m                 \u001b[0margs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mmagic_arg_s\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcell\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2371\u001b[0;31m                 \u001b[0mresult\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2372\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0mresult\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2373\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/conda/lib/python3.6/site-packages/IPython/core/magics/script.py\u001b[0m in \u001b[0;36mnamed_script_magic\u001b[0;34m(line, cell)\u001b[0m\n\u001b[1;32m    140\u001b[0m             \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    141\u001b[0m                 \u001b[0mline\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mscript\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 142\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshebang\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mline\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcell\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    143\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    144\u001b[0m         \u001b[0;31m# write a basic docstring:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<decorator-gen-110>\u001b[0m in \u001b[0;36mshebang\u001b[0;34m(self, line, cell)\u001b[0m\n",
      "\u001b[0;32m/opt/conda/lib/python3.6/site-packages/IPython/core/magic.py\u001b[0m in \u001b[0;36m<lambda>\u001b[0;34m(f, *a, **k)\u001b[0m\n\u001b[1;32m    185\u001b[0m     \u001b[0;31m# but it's overkill for just that one bit of state.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    186\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mmagic_deco\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 187\u001b[0;31m         \u001b[0mcall\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mlambda\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    188\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    189\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mcallable\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/conda/lib/python3.6/site-packages/IPython/core/magics/script.py\u001b[0m in \u001b[0;36mshebang\u001b[0;34m(self, line, cell)\u001b[0m\n\u001b[1;32m    243\u001b[0m             \u001b[0msys\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstderr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mflush\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    244\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mraise_error\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreturncode\u001b[0m\u001b[0;34m!=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 245\u001b[0;31m             \u001b[0;32mraise\u001b[0m \u001b[0mCalledProcessError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreturncode\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcell\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0moutput\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mout\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstderr\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0merr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    246\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    247\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m_run_script\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mp\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcell\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mto_close\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mCalledProcessError\u001b[0m: Command 'b'dx download dbVar_byAssembly_GRCh37/GRCh37.variant_call.all.vcf.gz\\n'' returned non-zero exit status 1."
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "dx download dbVar_byAssembly_GRCh37/GRCh37.variant_call.all.vcf.gz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "outf = open('all_variants_chr21_grch37.tsv', 'w')\n",
    "# header\n",
    "outf.write('chr\\tstart\\tend\\ttype\\tvariant_id\\n')\n",
    "# read VCF with all variants\n",
    "svtypes = {}\n",
    "for line in gzip.open('GRCh37.variant_call.all.vcf.gz', 'rb'):\n",
    "    line = line.decode('ascii')\n",
    "    # skip headers\n",
    "    if line[0] == '#':\n",
    "        continue\n",
    "    line = line.rstrip().split('\\t')\n",
    "    # skip if not chr 21\n",
    "    if line[0] != '21':\n",
    "        continue\n",
    "    # parse INFO field\n",
    "    infos = {}\n",
    "    for info in line[7].split(';'):\n",
    "        info_pair = info.split('=')\n",
    "        if len(info_pair) == 1:\n",
    "            infos[info] = True\n",
    "        else:\n",
    "            infos[info_pair[0]] = info_pair[1]\n",
    "    # prepare info for TSV\n",
    "    out_line = [line[0]]\n",
    "    # start pos\n",
    "    pos_s = line[1]\n",
    "    out_line.append(pos_s)\n",
    "    # end pos\n",
    "    pos_e = pos_s\n",
    "    if 'END' in infos:\n",
    "        pos_e = infos['END']\n",
    "    out_line.append(pos_e)\n",
    "    # SV type\n",
    "    svtype = infos['SVTYPE']\n",
    "    out_line.append(svtype)\n",
    "    svtypes[svtype] = True\n",
    "    # skip if BND type\n",
    "    if svtype == 'BND':\n",
    "        continue\n",
    "    if int(pos_e) < int(pos_s):\n",
    "        print(svtype)\n",
    "    # variant id\n",
    "    out_line.append(line[2])\n",
    "    # write line in tsv\n",
    "    outf.write('\\t'.join(out_line) + '\\n')\n",
    "outf.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dict_keys(['DEL', 'DUP', 'CNV', 'INS', 'BND', 'INV'])\n"
     ]
    }
   ],
   "source": [
    "print(svtypes.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ID                  file-Fy33J480Z5b09bvV4jpxj09z\n",
      "Class               file\n",
      "Project             project-Fy1b7V80Z5b4jXb224P1fY4b\n",
      "Folder              /chr-subset-genes-variants\n",
      "Name                all_variants_chr21_grch37.tsv\n",
      "State               closing\n",
      "Visibility          visible\n",
      "Types               -\n",
      "Properties          -\n",
      "Tags                -\n",
      "Outgoing links      -\n",
      "Created             Tue Oct 13 23:21:53 2020\n",
      "Created by          jmonlong\n",
      " via the job        job-Fy332xj0Z5bKZz0F8v9gKVfJ\n",
      "Last modified       Tue Oct 13 23:21:53 2020\n",
      "Media type          \n",
      "archivalState       \"live\"\n",
      "cloudAccount        \"cloudaccount-dnanexus\"\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "dx upload all_variants_chr21_grch37.tsv --path chr-subset-genes-variants/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extract genes in chr 21\n",
    "Gene names as gene IDs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "dx download Annotation/gene/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read VCF with all variants\n",
    "genes = {}\n",
    "for line in gzip.open('gencode.v19.annotation.gff3.gz', 'rb'):\n",
    "    line = line.decode('ascii')\n",
    "    # skip headers\n",
    "    if line[0] == '#':\n",
    "        continue\n",
    "    line = line.rstrip().split('\\t')\n",
    "    # skip if not chr 21 or the \"gene\" info\n",
    "    if line[0] != 'chr21' or line[2] != 'gene':\n",
    "        continue\n",
    "    # parse INFO field\n",
    "    infos = {}\n",
    "    for info in line[8].split(';'):\n",
    "        info_pair = info.split('=')\n",
    "        if len(info_pair) == 1:\n",
    "            infos[info] = True\n",
    "        else:\n",
    "            infos[info_pair[0]] = info_pair[1]\n",
    "    if infos['gene_type'] == 'protein_coding':\n",
    "        genes[infos['gene_name']] = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "241"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(genes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "outf = open('genes_chr21.tsv', 'w')\n",
    "# header\n",
    "outf.write('gene_id\\n')\n",
    "for gene in genes.keys():\n",
    "    outf.write(gene + '\\n')\n",
    "outf.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ID                  file-Fy31xB00Z5bFYXyG5Bf33yzf\n",
      "Class               file\n",
      "Project             project-Fy1b7V80Z5b4jXb224P1fY4b\n",
      "Folder              /chr-subset-genes-variants\n",
      "Name                genes_chr21.tsv\n",
      "State               closing\n",
      "Visibility          visible\n",
      "Types               -\n",
      "Properties          -\n",
      "Tags                -\n",
      "Outgoing links      -\n",
      "Created             Tue Oct 13 21:39:52 2020\n",
      "Created by          jmonlong\n",
      " via the job        job-Fy30f5Q0Z5b2Pk7kByqz8fk6\n",
      "Last modified       Tue Oct 13 21:39:52 2020\n",
      "Media type          \n",
      "archivalState       \"live\"\n",
      "cloudAccount        \"cloudaccount-dnanexus\"\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "dx upload genes_chr21.tsv --path chr-subset-genes-variants/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Save the notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ID                  file-Fy31xV80Z5bFYXyG5Bf33z02\n",
      "Class               file\n",
      "Project             project-Fy1b7V80Z5b4jXb224P1fY4b\n",
      "Folder              /chr-subset-genes-variants\n",
      "Name                extract-chr21-genes-variants.ipynb\n",
      "State               closing\n",
      "Visibility          visible\n",
      "Types               -\n",
      "Properties          -\n",
      "Tags                -\n",
      "Outgoing links      -\n",
      "Created             Tue Oct 13 21:40:21 2020\n",
      "Created by          jmonlong\n",
      " via the job        job-Fy30f5Q0Z5b2Pk7kByqz8fk6\n",
      "Last modified       Tue Oct 13 21:40:22 2020\n",
      "Media type          \n",
      "archivalState       \"live\"\n",
      "cloudAccount        \"cloudaccount-dnanexus\"\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "dx upload extract-chr21-genes-variants.ipynb --path chr-subset-genes-variants/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}