dysgu_api_demo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dysgu API for python3\n",
    "==================\n",
    "\n",
    "This notebook gives an overview of the dysgu-api for python.\n",
    "\n",
    "\n",
    "Calling structural variants\n",
    "------------------------------\n",
    "\n",
    "Dysgu relies on the `pysam` library to provide access to alignment files and the reference genome. To call SVs, first initialize a `DysguSV` class with an indexed reference genome and alignment file using pysam. Here we will use the HG002 samples which is ~40x coverage:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "import pysam\n",
    "from dysgu import DysguSV\n",
    "\n",
    "bam = pysam.AlignmentFile('HG002.bam', 'rb')\n",
    "genome = pysam.FastaFile('ucsc.hg19.fasta')\n",
    "\n",
    "dysgu = DysguSV(genome, bam)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calling SVs is performed by passing an iterable of reads to dysgu. For example to call SVs over the first 1 Mb of chr1:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time: 1.890 seconds\n"
     ]
    }
   ],
   "source": [
    "t0 = time.time()\n",
    "\n",
    "df = dysgu(\"chr1:1-1000000\")\n",
    "\n",
    "print(f\"Time: {time.time() - t0:.3f} seconds\")\n",
    "\n",
    "# Dysgu can also be called with an iterable, for example:\n",
    "f_iter = bam.fetch(\"chr1\", 10_000, 20_000)\n",
    "df2 = dysgu(f_iter)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dysgu will return a pandas DataFrame with the resulting calls, or `None` if no SVs were identified:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>chrA</th>\n",
       "      <th>posA</th>\n",
       "      <th>event_id</th>\n",
       "      <th>ref_seq</th>\n",
       "      <th>variant_seq</th>\n",
       "      <th>filter</th>\n",
       "      <th>sample</th>\n",
       "      <th>svtype</th>\n",
       "      <th>posB</th>\n",
       "      <th>chrB</th>\n",
       "      <th>...</th>\n",
       "      <th>inner_cn</th>\n",
       "      <th>outer_cn</th>\n",
       "      <th>compress</th>\n",
       "      <th>ref_rep</th>\n",
       "      <th>jitter</th>\n",
       "      <th>contigA</th>\n",
       "      <th>right_ins_seq</th>\n",
       "      <th>left_ins_seq</th>\n",
       "      <th>partners</th>\n",
       "      <th>prob</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>chr1</td>\n",
       "      <td>54720</td>\n",
       "      <td>76</td>\n",
       "      <td>C</td>\n",
       "      <td>TTTCTTTCTTTTTTTCTTTTTTTTTTTCTTTCTTTCTTTCTTTCTT...</td>\n",
       "      <td>PASS</td>\n",
       "      <td>sample</td>\n",
       "      <td>INS</td>\n",
       "      <td>54720</td>\n",
       "      <td>chr1</td>\n",
       "      <td>...</td>\n",
       "      <td>2.25</td>\n",
       "      <td>0.000</td>\n",
       "      <td>31.219999</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>ccttctattttttttttctttctttctttttttctttttttttttc...</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>0.626</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>chr1</td>\n",
       "      <td>66321</td>\n",
       "      <td>81</td>\n",
       "      <td>A</td>\n",
       "      <td>&lt;INS&gt;</td>\n",
       "      <td>lowProb</td>\n",
       "      <td>sample</td>\n",
       "      <td>INS</td>\n",
       "      <td>66321</td>\n",
       "      <td>chr1</td>\n",
       "      <td>...</td>\n",
       "      <td>16.50</td>\n",
       "      <td>1.500</td>\n",
       "      <td>19.750000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>aatatatatattattatataatataatatatattATATAATATATA...</td>\n",
       "      <td>None</td>\n",
       "      <td>aatatatatattattatataatataatatatatt</td>\n",
       "      <td>None</td>\n",
       "      <td>0.427</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>chr1</td>\n",
       "      <td>66452</td>\n",
       "      <td>85</td>\n",
       "      <td>A</td>\n",
       "      <td>&lt;INS&gt;</td>\n",
       "      <td>PASS</td>\n",
       "      <td>sample</td>\n",
       "      <td>INS</td>\n",
       "      <td>66452</td>\n",
       "      <td>chr1</td>\n",
       "      <td>...</td>\n",
       "      <td>17.50</td>\n",
       "      <td>0.750</td>\n",
       "      <td>19.020000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>ATATTATTATATAATATAATATATATTATATAATATATTTTATTAT...</td>\n",
       "      <td>atatatattatataatatattttattatataaatatatattatattat</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>0.758</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>chr1</td>\n",
       "      <td>724219</td>\n",
       "      <td>187</td>\n",
       "      <td>A</td>\n",
       "      <td>TGGAATGAAGTGGAATGGAGTGGAATTGATTGGTGCGGAGCAGAAA...</td>\n",
       "      <td>lowProb</td>\n",
       "      <td>sample</td>\n",
       "      <td>INS</td>\n",
       "      <td>724219</td>\n",
       "      <td>chr1</td>\n",
       "      <td>...</td>\n",
       "      <td>52.25</td>\n",
       "      <td>21.750</td>\n",
       "      <td>36.560001</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>tggaatgcaatgaaacagaatggaatggattggatttgtatggact...</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>0.339</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>chr1</td>\n",
       "      <td>724350</td>\n",
       "      <td>199</td>\n",
       "      <td>T</td>\n",
       "      <td>&lt;DEL&gt;</td>\n",
       "      <td>lowProb</td>\n",
       "      <td>sample</td>\n",
       "      <td>DEL</td>\n",
       "      <td>724468</td>\n",
       "      <td>chr1</td>\n",
       "      <td>...</td>\n",
       "      <td>7.00</td>\n",
       "      <td>12.375</td>\n",
       "      <td>36.139999</td>\n",
       "      <td>0.320301</td>\n",
       "      <td>0.0</td>\n",
       "      <td></td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>0.076</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>chr1</td>\n",
       "      <td>724605</td>\n",
       "      <td>231</td>\n",
       "      <td>T</td>\n",
       "      <td>AAGGAATTTATTCGAATGGAATAGAATGGAAAGGA</td>\n",
       "      <td>PASS</td>\n",
       "      <td>sample</td>\n",
       "      <td>INS</td>\n",
       "      <td>724708</td>\n",
       "      <td>chr1</td>\n",
       "      <td>...</td>\n",
       "      <td>53.50</td>\n",
       "      <td>22.000</td>\n",
       "      <td>37.990002</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>ccGAATGGAATGGAATGTAATGGAACGGAATGGAATGGAACAGAAT...</td>\n",
       "      <td>tgaatggaaatgaaagtaata</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>0.518</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>chr1</td>\n",
       "      <td>725033</td>\n",
       "      <td>196</td>\n",
       "      <td>T</td>\n",
       "      <td>&lt;INS&gt;</td>\n",
       "      <td>lowProb</td>\n",
       "      <td>sample</td>\n",
       "      <td>INS</td>\n",
       "      <td>726715</td>\n",
       "      <td>chr1</td>\n",
       "      <td>...</td>\n",
       "      <td>22.50</td>\n",
       "      <td>7.000</td>\n",
       "      <td>37.020000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td></td>\n",
       "      <td>None</td>\n",
       "      <td>atgtacttgaatggagtagaatggaatggaatgtaaaggaaatgaa...</td>\n",
       "      <td>None</td>\n",
       "      <td>0.116</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>chr1</td>\n",
       "      <td>725302</td>\n",
       "      <td>318</td>\n",
       "      <td>G</td>\n",
       "      <td>&lt;INS&gt;</td>\n",
       "      <td>lowProb</td>\n",
       "      <td>sample</td>\n",
       "      <td>INS</td>\n",
       "      <td>725361</td>\n",
       "      <td>chr1</td>\n",
       "      <td>...</td>\n",
       "      <td>20.25</td>\n",
       "      <td>21.750</td>\n",
       "      <td>33.680000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>tgaaatggaatggaatggaatagattccaaagaatggaatggactc...</td>\n",
       "      <td>None</td>\n",
       "      <td>tgaaatggaatggaatggaatagattccaaagaatggaatggactc...</td>\n",
       "      <td>None</td>\n",
       "      <td>0.286</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>chr1</td>\n",
       "      <td>725717</td>\n",
       "      <td>351</td>\n",
       "      <td>G</td>\n",
       "      <td>&lt;INS&gt;</td>\n",
       "      <td>lowProb</td>\n",
       "      <td>sample</td>\n",
       "      <td>INS</td>\n",
       "      <td>725763</td>\n",
       "      <td>chr1</td>\n",
       "      <td>...</td>\n",
       "      <td>46.00</td>\n",
       "      <td>13.500</td>\n",
       "      <td>35.900002</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>ggaagggaatggagtgtacaagaatggaatgggatggcaaaaaatg...</td>\n",
       "      <td>None</td>\n",
       "      <td>ggaagggaatggagtgtacaagaatggaatgggatggcaaaaaatg...</td>\n",
       "      <td>None</td>\n",
       "      <td>0.211</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>chr1</td>\n",
       "      <td>726095</td>\n",
       "      <td>358</td>\n",
       "      <td>T</td>\n",
       "      <td>&lt;DEL&gt;</td>\n",
       "      <td>lowProb</td>\n",
       "      <td>sample</td>\n",
       "      <td>DEL</td>\n",
       "      <td>726125</td>\n",
       "      <td>chr1</td>\n",
       "      <td>...</td>\n",
       "      <td>10.00</td>\n",
       "      <td>7.875</td>\n",
       "      <td>41.450001</td>\n",
       "      <td>0.305509</td>\n",
       "      <td>0.0</td>\n",
       "      <td>GGAATAGAATAGAATGGACTGGAATGTAATGAGTTTGGAATGGAAT...</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>0.121</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>10 rows × 78 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   chrA    posA  event_id ref_seq  \\\n",
       "0  chr1   54720        76       C   \n",
       "1  chr1   66321        81       A   \n",
       "2  chr1   66452        85       A   \n",
       "3  chr1  724219       187       A   \n",
       "4  chr1  724350       199       T   \n",
       "5  chr1  724605       231       T   \n",
       "6  chr1  725033       196       T   \n",
       "7  chr1  725302       318       G   \n",
       "8  chr1  725717       351       G   \n",
       "9  chr1  726095       358       T   \n",
       "\n",
       "                                         variant_seq   filter  sample svtype  \\\n",
       "0  TTTCTTTCTTTTTTTCTTTTTTTTTTTCTTTCTTTCTTTCTTTCTT...     PASS  sample    INS   \n",
       "1                                              <INS>  lowProb  sample    INS   \n",
       "2                                              <INS>     PASS  sample    INS   \n",
       "3  TGGAATGAAGTGGAATGGAGTGGAATTGATTGGTGCGGAGCAGAAA...  lowProb  sample    INS   \n",
       "4                                              <DEL>  lowProb  sample    DEL   \n",
       "5                AAGGAATTTATTCGAATGGAATAGAATGGAAAGGA     PASS  sample    INS   \n",
       "6                                              <INS>  lowProb  sample    INS   \n",
       "7                                              <INS>  lowProb  sample    INS   \n",
       "8                                              <INS>  lowProb  sample    INS   \n",
       "9                                              <DEL>  lowProb  sample    DEL   \n",
       "\n",
       "     posB  chrB  ...  inner_cn  outer_cn   compress   ref_rep  jitter  \\\n",
       "0   54720  chr1  ...      2.25     0.000  31.219999  0.000000     0.0   \n",
       "1   66321  chr1  ...     16.50     1.500  19.750000  0.000000     0.0   \n",
       "2   66452  chr1  ...     17.50     0.750  19.020000  0.000000     0.0   \n",
       "3  724219  chr1  ...     52.25    21.750  36.560001  0.000000     0.0   \n",
       "4  724468  chr1  ...      7.00    12.375  36.139999  0.320301     0.0   \n",
       "5  724708  chr1  ...     53.50    22.000  37.990002  0.000000     0.0   \n",
       "6  726715  chr1  ...     22.50     7.000  37.020000  0.000000     0.0   \n",
       "7  725361  chr1  ...     20.25    21.750  33.680000  0.000000     0.0   \n",
       "8  725763  chr1  ...     46.00    13.500  35.900002  0.000000     0.0   \n",
       "9  726125  chr1  ...     10.00     7.875  41.450001  0.305509     0.0   \n",
       "\n",
       "                                             contigA  \\\n",
       "0  ccttctattttttttttctttctttctttttttctttttttttttc...   \n",
       "1  aatatatatattattatataatataatatatattATATAATATATA...   \n",
       "2  ATATTATTATATAATATAATATATATTATATAATATATTTTATTAT...   \n",
       "3  tggaatgcaatgaaacagaatggaatggattggatttgtatggact...   \n",
       "4                                                      \n",
       "5  ccGAATGGAATGGAATGTAATGGAACGGAATGGAATGGAACAGAAT...   \n",
       "6                                                      \n",
       "7  tgaaatggaatggaatggaatagattccaaagaatggaatggactc...   \n",
       "8  ggaagggaatggagtgtacaagaatggaatgggatggcaaaaaatg...   \n",
       "9  GGAATAGAATAGAATGGACTGGAATGTAATGAGTTTGGAATGGAAT...   \n",
       "\n",
       "                                      right_ins_seq  \\\n",
       "0                                              None   \n",
       "1                                              None   \n",
       "2  atatatattatataatatattttattatataaatatatattatattat   \n",
       "3                                              None   \n",
       "4                                              None   \n",
       "5                             tgaatggaaatgaaagtaata   \n",
       "6                                              None   \n",
       "7                                              None   \n",
       "8                                              None   \n",
       "9                                              None   \n",
       "\n",
       "                                        left_ins_seq  partners   prob  \n",
       "0                                               None      None  0.626  \n",
       "1                 aatatatatattattatataatataatatatatt      None  0.427  \n",
       "2                                               None      None  0.758  \n",
       "3                                               None      None  0.339  \n",
       "4                                               None      None  0.076  \n",
       "5                                               None      None  0.518  \n",
       "6  atgtacttgaatggagtagaatggaatggaatgtaaaggaaatgaa...      None  0.116  \n",
       "7  tgaaatggaatggaatggaatagattccaaagaatggaatggactc...      None  0.286  \n",
       "8  ggaagggaatggagtgtacaagaatggaatgggatggcaaaaaatg...      None  0.211  \n",
       "9                                               None      None  0.121  \n",
       "\n",
       "[10 rows x 78 columns]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.head(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "lowProb    15\n",
      "PASS       12\n",
      "Name: filter, dtype: int64\n",
      "INS    9\n",
      "DEL    3\n",
      "Name: svtype, dtype: int64\n"
     ]
    }
   ],
   "source": [
    "# display the number of calls with a PASS in the filter column:\n",
    "print(df['filter'].value_counts())\n",
    "\n",
    "# display the counts of SV types with a PASS\n",
    "print(df.query(\" filter == 'PASS' \")['svtype'].value_counts());"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notes\n",
    "-------\n",
    "- Dysgu expects the input alignment file to be in sorted order\n",
    "\n",
    "- There is currently some overhead associated with preparing a new chromosome for SV calling, for this reason it will be inefficient to repeatedly call multiple small regions. It is far more efficient to iterate over all target regions in one go. Dysgu provides a few tools for manipulating bed files to support this pattern:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[('chr1', 10000, 21000), ('chr1', 25000, 30000), ('chr1', 35000, 36000), ('chr1', 1500000, 16000000)]\n",
      "11.701174020767212\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "lowProb    189\n",
       "PASS       132\n",
       "Name: filter, dtype: int64"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from dysgu import load_bed, merge_intervals\n",
    "\n",
    "bed = load_bed('example.bed')\n",
    "\n",
    "# merge and sort intervals\n",
    "bed = merge_intervals(bed, srt=True)\n",
    "print(bed)\n",
    "\n",
    "t0 = time.time()\n",
    "\n",
    "df = dysgu.call_bed_regions(bed)\n",
    "\n",
    "print(time.time() - t0)\n",
    "df['filter'].value_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Options\n",
    "---------\n",
    "\n",
    "The default options for dysgu are similar to the CLI application. These can be inspected using:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'add_gt': True,\n",
      " 'add_kind': True,\n",
      " 'all_sites': False,\n",
      " 'buffer_size': 10000,\n",
      " 'clean': False,\n",
      " 'clip_length': 15,\n",
      " 'contigs': True,\n",
      " 'diploid': True,\n",
      " 'dist_norm': 100,\n",
      " 'drop_gaps': True,\n",
      " 'ibam': None,\n",
      " 'keep_small': False,\n",
      " 'low_mem': False,\n",
      " 'max_cov': 200,\n",
      " 'max_tlen': 1000,\n",
      " 'merge_dist': None,\n",
      " 'merge_within': True,\n",
      " 'metrics': True,\n",
      " 'min_size': 30,\n",
      " 'min_support': 3,\n",
      " 'mode': 'pe',\n",
      " 'model': None,\n",
      " 'mq': 1,\n",
      " 'overwrite': True,\n",
      " 'paired': True,\n",
      " 'parse_probs': False,\n",
      " 'pfix': 'dysgu_reads',\n",
      " 'pl': 'pe',\n",
      " 'procs': 1,\n",
      " 'reference': None,\n",
      " 'regions': None,\n",
      " 'regions_mm_only': False,\n",
      " 'regions_only': False,\n",
      " 'remap': True,\n",
      " 'sites': None,\n",
      " 'sites_pass_only': True,\n",
      " 'sites_prob': 0.6,\n",
      " 'spd': 0.3,\n",
      " 'sv_aligns': None,\n",
      " 'template_size': '',\n",
      " 'thresholds': {'DEL': 0.45,\n",
      "                'DUP': 0.45,\n",
      "                'INS': 0.45,\n",
      "                'INV': 0.45,\n",
      "                'TRA': 0.45},\n",
      " 'trust_ins_len': True,\n",
      " 'verbosity': 2,\n",
      " 'working_directory': 'tempfile',\n",
      " 'z_breadth': 2,\n",
      " 'z_depth': 2}\n"
     ]
    }
   ],
   "source": [
    "import pprint\n",
    "from dysgu import dysgu_default_args\n",
    "\n",
    "pprint.pprint(dysgu_default_args())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Preset options are also available using the `dysgu_preset_args(mode)` function, where mode can be \"pacbio-revio\" etc.\n",
    "\n",
    "Options can be specified as key word arguments during initialization. For example, to change the `min-support` option:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "PASS       87\n",
       "lowProb    87\n",
       "Name: filter, dtype: int64"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dysgu = DysguSV(genome, bam, min_support=5)\n",
    "df = dysgu.call_bed_regions(bed)\n",
    "df['filter'].value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "PASS       46\n",
       "lowProb    34\n",
       "Name: filter, dtype: int64"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Or alternatively:\n",
    "dysgu.set_option(\"min_support\", 10)\n",
    "\n",
    "dysgu.call_bed_regions(bed)['filter'].value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "PASS       26\n",
       "lowProb    14\n",
       "Name: filter, dtype: int64"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Or provide a mapping of arguments:\n",
    "dysgu.set_option({\"min_support\": 10, \"mq\": 20, \"min_size\": 100})\n",
    "\n",
    "dysgu.call_bed_regions(bed)['filter'].value_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Saving and loading a VCF\n",
    "------------------------------\n",
    "\n",
    "`to_vcf` and `load_vcf` functions can be used as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dysgu import load_dysgu_vcf\n",
    "import pandas as pd\n",
    "\n",
    "\n",
    "dysgu = DysguSV(genome, bam)\n",
    "\n",
    "df = dysgu(\"chr1:1000000-2000000\")\n",
    "\n",
    "# Note a sample name in needed for the vcf file, this can be set as an attribute,\n",
    "# or during construction i.e. DysguSV(genome, bam, 'sample1')\n",
    "dysgu.sample_name = \"sample1\"\n",
    "\n",
    "with open(\"output.vcf\", \"w\") as out:\n",
    "    dysgu.to_vcf(df, out)\n",
    "\n",
    "df2 = load_dysgu_vcf(\"output.vcf\")\n",
    "\n",
    "\n",
    "# just to check dataframes are the same. Note float values are rounded when saving to vcf\n",
    "pd.testing.assert_frame_equal(df[[\"chrA\", \"posA\", \"event_id\"]],\n",
    "                              df2[[\"chrA\", \"posA\", \"event_id\"]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Switching machine learning models\n",
    "------------------------------------------\n",
    "\n",
    "To change the ML model, first configure dysgu with the desired options, then use the `apply_model` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "lowProb    21\n",
      "PASS       16\n",
      "Name: filter, dtype: int64\n",
      "lowProb    24\n",
      "PASS       13\n",
      "Name: filter, dtype: int64\n",
      "lowProb    23\n",
      "PASS       14\n",
      "Name: filter, dtype: int64\n",
      "PASS       19\n",
      "lowProb    18\n",
      "Name: filter, dtype: int64\n"
     ]
    }
   ],
   "source": [
    "dysgu = DysguSV(genome, bam)\n",
    "df = dysgu(\"chr1:1000000-2000000\")\n",
    "\n",
    "# try the non-diploid model\n",
    "dysgu.set_option(\"diploid\", False)\n",
    "df_non_diploid = dysgu.apply_model(df)\n",
    "print(df_non_diploid['filter'].value_counts())\n",
    "\n",
    "# try the diploid model\n",
    "dysgu.set_option(\"diploid\", True)\n",
    "df_diploid = dysgu.apply_model(df)\n",
    "print(df_diploid['filter'].value_counts())\n",
    "\n",
    "# try the no-contigs model\n",
    "dysgu.set_option(\"contigs\", False)\n",
    "df_no_conts = dysgu.apply_model(df)\n",
    "print(df_no_conts['filter'].value_counts())\n",
    "\n",
    "# try the pacbio model (why not)\n",
    "dysgu.set_option({\"contigs\": False, \"pl\": \"pacbio\"})\n",
    "df_no_conts = dysgu.apply_model(df)\n",
    "print(df_no_conts['filter'].value_counts())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Merging SVs\n",
    "---------------\n",
    "\n",
    "Merging is supported via the `merge_dysgu_df` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dysgu import merge_dysgu_df\n",
    "\n",
    "\n",
    "samp1 = DysguSV(genome, pysam.AlignmentFile('HG002.bam', 'rb'), 'HG002')\n",
    "samp2 = DysguSV(genome, pysam.AlignmentFile('NA12878.bwa.bam', 'rb'), 'NA12878')\n",
    "\n",
    "r = \"chr1:8705980-8706344\"\n",
    "\n",
    "df_samp1 = samp1(r)\n",
    "df_samp2 = samp2(r)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>chrA</th>\n",
       "      <th>posA</th>\n",
       "      <th>event_id</th>\n",
       "      <th>ref_seq</th>\n",
       "      <th>variant_seq</th>\n",
       "      <th>filter</th>\n",
       "      <th>sample</th>\n",
       "      <th>svtype</th>\n",
       "      <th>posB</th>\n",
       "      <th>chrB</th>\n",
       "      <th>...</th>\n",
       "      <th>inner_cn</th>\n",
       "      <th>outer_cn</th>\n",
       "      <th>compress</th>\n",
       "      <th>ref_rep</th>\n",
       "      <th>jitter</th>\n",
       "      <th>contigA</th>\n",
       "      <th>right_ins_seq</th>\n",
       "      <th>left_ins_seq</th>\n",
       "      <th>partners</th>\n",
       "      <th>prob</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>chr1</td>\n",
       "      <td>8706162</td>\n",
       "      <td>1</td>\n",
       "      <td>T</td>\n",
       "      <td>&lt;DEL&gt;</td>\n",
       "      <td>PASS</td>\n",
       "      <td>HG002</td>\n",
       "      <td>DEL</td>\n",
       "      <td>8706208</td>\n",
       "      <td>chr1</td>\n",
       "      <td>...</td>\n",
       "      <td>1.7</td>\n",
       "      <td>0.0</td>\n",
       "      <td>47.43</td>\n",
       "      <td>0.183956</td>\n",
       "      <td>0.0</td>\n",
       "      <td>ACTATAATGTGACTGAAAATAAGAACTACATTTTCCCAACTGTTCC...</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>0.887</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1 rows × 78 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   chrA     posA  event_id ref_seq variant_seq filter sample svtype     posB  \\\n",
       "0  chr1  8706162         1       T       <DEL>   PASS  HG002    DEL  8706208   \n",
       "\n",
       "   chrB  ...  inner_cn  outer_cn compress   ref_rep  jitter  \\\n",
       "0  chr1  ...       1.7       0.0    47.43  0.183956     0.0   \n",
       "\n",
       "                                             contigA right_ins_seq  \\\n",
       "0  ACTATAATGTGACTGAAAATAAGAACTACATTTTCCCAACTGTTCC...          None   \n",
       "\n",
       "  left_ins_seq  partners   prob  \n",
       "0         None      None  0.887  \n",
       "\n",
       "[1 rows x 78 columns]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_samp1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>chrA</th>\n",
       "      <th>posA</th>\n",
       "      <th>event_id</th>\n",
       "      <th>ref_seq</th>\n",
       "      <th>variant_seq</th>\n",
       "      <th>filter</th>\n",
       "      <th>sample</th>\n",
       "      <th>svtype</th>\n",
       "      <th>posB</th>\n",
       "      <th>chrB</th>\n",
       "      <th>...</th>\n",
       "      <th>inner_cn</th>\n",
       "      <th>outer_cn</th>\n",
       "      <th>compress</th>\n",
       "      <th>ref_rep</th>\n",
       "      <th>jitter</th>\n",
       "      <th>contigA</th>\n",
       "      <th>right_ins_seq</th>\n",
       "      <th>left_ins_seq</th>\n",
       "      <th>partners</th>\n",
       "      <th>prob</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>chr1</td>\n",
       "      <td>8706162</td>\n",
       "      <td>1</td>\n",
       "      <td>T</td>\n",
       "      <td>&lt;DEL&gt;</td>\n",
       "      <td>PASS</td>\n",
       "      <td>NA12878</td>\n",
       "      <td>DEL</td>\n",
       "      <td>8706208</td>\n",
       "      <td>chr1</td>\n",
       "      <td>...</td>\n",
       "      <td>0.833333</td>\n",
       "      <td>0.0</td>\n",
       "      <td>47.939999</td>\n",
       "      <td>0.183956</td>\n",
       "      <td>0.0</td>\n",
       "      <td>CACTATAATGTGACTGAAAATAAGAACTACATTTTCCCAACTGTTC...</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>0.795</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1 rows × 78 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   chrA     posA  event_id ref_seq variant_seq filter   sample svtype  \\\n",
       "0  chr1  8706162         1       T       <DEL>   PASS  NA12878    DEL   \n",
       "\n",
       "      posB  chrB  ...  inner_cn  outer_cn   compress   ref_rep  jitter  \\\n",
       "0  8706208  chr1  ...  0.833333       0.0  47.939999  0.183956     0.0   \n",
       "\n",
       "                                             contigA right_ins_seq  \\\n",
       "0  CACTATAATGTGACTGAAAATAAGAACTACATTTTCCCAACTGTTC...          None   \n",
       "\n",
       "  left_ins_seq  partners   prob  \n",
       "0         None      None  0.795  \n",
       "\n",
       "[1 rows x 78 columns]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_samp2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>chrA</th>\n",
       "      <th>posA</th>\n",
       "      <th>event_id</th>\n",
       "      <th>ref_seq</th>\n",
       "      <th>variant_seq</th>\n",
       "      <th>filter</th>\n",
       "      <th>sample</th>\n",
       "      <th>svtype</th>\n",
       "      <th>posB</th>\n",
       "      <th>chrB</th>\n",
       "      <th>...</th>\n",
       "      <th>outer_cn</th>\n",
       "      <th>compress</th>\n",
       "      <th>ref_rep</th>\n",
       "      <th>jitter</th>\n",
       "      <th>contigA</th>\n",
       "      <th>contigA</th>\n",
       "      <th>right_ins_seq</th>\n",
       "      <th>left_ins_seq</th>\n",
       "      <th>partners</th>\n",
       "      <th>prob</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>chr1</td>\n",
       "      <td>8706162</td>\n",
       "      <td>0</td>\n",
       "      <td>T</td>\n",
       "      <td>&lt;DEL&gt;</td>\n",
       "      <td>PASS</td>\n",
       "      <td>HG002</td>\n",
       "      <td>DEL</td>\n",
       "      <td>8706208</td>\n",
       "      <td>chr1</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>47.43</td>\n",
       "      <td>0.183956</td>\n",
       "      <td>0.0</td>\n",
       "      <td>ACTATAATGTGACTGAAAATAAGAACTACATTTTCCCAACTGTTCC...</td>\n",
       "      <td>ACTATAATGTGACTGAAAATAAGAACTACATTTTCCCAACTGTTCC...</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>[(1, 1)]</td>\n",
       "      <td>0.887</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1 rows × 80 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   chrA     posA  event_id ref_seq variant_seq filter sample svtype     posB  \\\n",
       "0  chr1  8706162         0       T       <DEL>   PASS  HG002    DEL  8706208   \n",
       "\n",
       "   chrB  ...  outer_cn  compress   ref_rep  jitter  \\\n",
       "0  chr1  ...       0.0     47.43  0.183956     0.0   \n",
       "\n",
       "                                             contigA  \\\n",
       "0  ACTATAATGTGACTGAAAATAAGAACTACATTTTCCCAACTGTTCC...   \n",
       "\n",
       "                                             contigA right_ins_seq  \\\n",
       "0  ACTATAATGTGACTGAAAATAAGAACTACATTTTCCCAACTGTTCC...          None   \n",
       "\n",
       "  left_ins_seq  partners   prob  \n",
       "0         None  [(1, 1)]  0.887  \n",
       "\n",
       "[1 rows x 80 columns]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "merged = merge_dysgu_df(df_samp1, df_samp2)\n",
    "merged"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Following merging, the `partners` column will contain a list of the merged variants. Each item merged is represented as a tuple with the format `(df_index, event_id)`. In the example above df_samp1 would have `df_index = 0` and df_samp2 would have `df_index = 1`. \n",
    "\n",
    "In the partners column, the list `[(1, 1)]` indicates that the variant was merged with a variant from df_index=1 with event_id=1, hence (1, 1).\n",
    "\n",
    "Note, the `sample` column is not modified during merging, so this is left to the user to manipulate if needs be"
   ]
  }
 ],
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}