HW3.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "little-appeal",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pysam\n",
    "from Bio import SeqIO\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from collections import Counter\n",
    "from tqdm import tqdm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "mental-tournament",
   "metadata": {},
   "source": [
    "# PacBio"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "european-clinic",
   "metadata": {},
   "source": [
    "##### >minimap2 -x map-pb -a -t 6 -o pacbio.sam reference.fasta pacbio_40x.fq.gz\n",
    "\n",
    "##### >samtools flagstat -@ 6 pacbio.sam > pacbio.stat\n",
    "\n",
    "52081 + 0 in total (QC-passed reads + QC-failed reads)\n",
    "\n",
    "1736 + 0 secondary\n",
    "\n",
    "8898 + 0 supplementary\n",
    "\n",
    "0 + 0 duplicates\n",
    "\n",
    "51199 + 0 mapped (98.31% : N/A)\n",
    "\n",
    "0 + 0 paired in sequencing\n",
    "\n",
    "0 + 0 read1\n",
    "\n",
    "0 + 0 read2\n",
    "\n",
    "0 + 0 properly paired (N/A : N/A)\n",
    "\n",
    "0 + 0 with itself and mate mapped\n",
    "\n",
    "0 + 0 singletons (N/A : N/A)\n",
    "\n",
    "0 + 0 with mate mapped to a different chr\n",
    "\n",
    "0 + 0 with mate mapped to a different chr (mapQ>=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "purple-roommate",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "========== pacbio_40x.fq.gz  ==========\n",
      "Number of substitutions:\n",
      "        A       C       G       T       -\n",
      "A       0   38461   30826   29043  303520\n",
      "C   29419       0   35823   26569  365946\n",
      "G   26566   35984       0   29173  397228\n",
      "T   29047   30698   37640       0  280733\n",
      "-  501939  505815  701449  515264       0\n",
      "mean_read_error=3.09937%\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2cAAAE/CAYAAADCCbvWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAfeUlEQVR4nO3dfZRtZ10f8O+PhDclSlKuWZckcIFGNCoN9Aq4sFahYqAuwZa2RMUUoelSWRWL1aBrKX2xoktEW60QCoLKiy/IkgLKu6W2CiYSIBBjAgkmMSShiNBqVcjTP86+5GQyM3dmzj5z9svns9asOWefM7Of57ef/ezznbPPnmqtBQAAgM2626YbAAAAgHAGAAAwCMIZAADAAAhnAAAAAyCcAQAADIBwBgAAMADCGQAAwAAIZwAcSFVdX1V/XVX327L8vVXVqurYhpoGAKMknAGwiuuSXHjiTlV9RZLP21xzNqeqTt10GwAYN+EMgFX8UpJvX7p/UZJfXH5CVd2zqn6yqv6kqm6pqhdV1b27x06vqjdU1W1V9Wfd7bOXfvZ3qurfV9X/rKpPV9Vbtr5Tt2VdT6qqK6rqU1X14aq6oFt+/6p6fVV9oqqurap/sbT8L6vqjKXf8fCq+nhV3b27/x1VdVXXvjdX1QOXntuq6rur6pok13TLfqaqbujacHlV/b2l59+7ql7R/a6rqur7q+rGpcfvX1Wv7epxXVX9q/1tDgDGTDgDYBW/n+QLqupLq+qUJE9N8stbnvP8JF+c5PwkfzvJWUl+uHvsbkl+IckDkzwgyV8m+dktP/8tSZ6e5IuS3CPJ923XkKp6ZBbB8N8kuW+Sr0lyfffwa5LcmOT+SZ6S5D9W1WNba3+a5PeS/OMt6/v11trfVNWTkvxgkn+U5EiS/5Hk1VtW/eQkj0pyXnf/D7q+npHkVUl+raru1T32I0mOJXlwkq9P8m1L7b9bkv+W5H1djR6X5NlV9Q3b9ReA6anW2qbbAMAIVdX1SZ6Z5NFJPj/Jf0/ynCRPSPI3SR6U5KNJ/k+Sh7XWPtz93FcleVVr7UHb/M7zk7yztXZ6d/93kryttfYfuvvfleSbWmsXbPOzL07yF621792y/JwsQtp9W2uf7pb9WJKjrbV/XlXPTPItrbXHVlUl+ZMk39pae1dV/VYWQe2l3c/drevPl7bWPlpVLcnjWmvv2KVOf5bka1tr76uqjyT5ztbam7vHnpnkea21s6vqUUl+rbX2gKWffW6SL26tPX2n3w/AdDg/HoBV/VKSd2URxn5xy2NHsvgM2uWL3JMkqSSnJElVfV6SFya5IMnp3eOnVdUprbXPdvc/tvT7/iLJfXZoxzlJ3rTN8vsn+cSJYNb5aJLj3e3XJvnPVXU0i3f4bs/iHbJk8Y7ez1TVC5Z+trJ4Z+uj3f0blldWVd+X5BndeluSL0hy4lTM+295/vLtBya5f1V9cmnZKUttAWDihDMAVtK9g3RdkidmEUqWfTyLUxW/rLV20zY//pwkD03yqNbax7p3zt6bRQDarxuSPGSb5X+a5IyqOm0poD0gyU1d+/+sqt6S5J8l+dIkr2l3nFZyQ5Ifba29cpf1fu4UlO7zZd+fxSmJH2yt3d69c3aiPzcnOTvJh7r752xp/3WttXP31FsAJsdnzgDowzOSPLa19n+XF7bWbk/ykiQvrKovSpKqOmvpc1SnZRHePtldlONHVmjDS5M8vaoeV1V369bzJa21G5L8ryQ/VlX3qqqHde1d/mzcq7K4sMlTutsnvCjJc6vqy7q2f2FV/ZNd2nBaks8kuS3JqVX1w1m8c3bCr3a/7/SqOivJs5Yee0+ST1fVD3QXDjmlqr68qr7yALUAYISEMwBW1lr7cGvtsh0e/oEk1yb5/ar6VJK3ZfFuWZL8dJJ7Z/EO2+8n+e0V2vCeLC4c8sIkf57FZ+BOXFnxwiwuxPGnSV6X5Edaa29b+vHXJzk3ycdaa+9b+p2vS/LjSV7Ttf3KLD5Tt5M3d3344yxOe/x/ufOpi/8uiwuTXJdFHX49yV916/pskm/M4mIi12VRk/+a5Av3XAQARs0FQQBgQ6rqO5M8tbX29zfdFgA2zztnAHBIqupoVT2mO+3yoVl85u51m24XAMPggiAAcHjukeTFWVzZ8pNZ/P+1/7LJBgEwHE5rBAAAGACnNQIAAAyAcAYAADAAh/qZs/vd737t2LFjh7lKAACAwbj88ss/3lo7st1jhxrOjh07lssu2+nf4AAAAExbVX10p8ec1ggAADAAwhkAAMAAnDScVdU5VfXOqvpQVX2wqr6nW/68qrqpqq7ovp64/uYCAABM014+c/aZJM9prf1hVZ2W5PKqemv32Atbaz+5vuYBAADMw0nDWWvt5iQ3d7c/XVVXJTlr3Q0DAACYk3195qyqjiV5eJJ3d4ueVVXvr6qXVdXpfTcOAABgLvYczqrqPklem+TZrbVPJfn5JA9Jcn4W76y9YIefu7iqLquqy2677bbVWwwAADBBewpnVXX3LILZK1trv5EkrbVbWmufba3dnuQlSR653c+21i5trR1vrR0/cmTb/7UGAAAwe3u5WmMleWmSq1prP7W0/OjS0745yZX9Nw8AAGAe9nK1xsckeVqSD1TVFd2yH0xyYVWdn6QluT7Jv1xD+wAAAGZhL1dr/N0ktc1Db+q/OQAAAPO0r6s1AgAAsB7CGQAAwAAIZwAAAAMgnAEAAAyAcAYAADAAwhkAAMAACGcAAAADIJwBAAAMgHAGAAAwAMIZAADAAAhnAAAAAyCcAQAADIBwBgAAMADCGQAAwAAIZwAAAAMgnAEAAAyAcAYAADAAwhkAAMAACGcAAAADIJwBAAAMgHAGAAAwAMIZAADAAAhnAAAAAyCcAQAADIBwBgAAMADCGQAAwAAIZwAAAAMgnAEAAAyAcAYAADAAwhkAAMAACGcAAAADIJwBAAAMgHAGAAAwAMIZAADAAAhnAAAAAyCcAQAADIBwBgAAMADCGQAAwAAIZwAAAAMgnAEAAAyAcAYAADAAwhkAAMAACGcAAGty7JI3broJwIgIZwAAAANw0nBWVedU1Tur6kNV9cGq+p5u+RlV9daquqb7fvr6mwsAADBNe3nn7DNJntNaOy/Jo5N8d1Wdl+SSJG9vrZ2b5O3dfQD4HKd0AX0xnzAHJw1nrbWbW2t/2N3+dJKrkpyV5ElJXtE97RVJnrymNgIAAEzevj5zVlXHkjw8ybuTnNlau7l76GNJztzhZy6uqsuq6rLbbrttlbYCAABM1p7DWVXdJ8lrkzy7tfap5cdaay1J2+7nWmuXttaOt9aOHzlyZKXGAgAATNWewllV3T2LYPbK1tpvdItvqaqj3eNHk9y6niYCAABM316u1lhJXprkqtbaTy099PokF3W3L0rym/03DwAAYB5O3cNzHpPkaUk+UFVXdMt+MMnzk/xqVT0jyUeT/NO1tBAAAGAGThrOWmu/m6R2ePhx/TYHAABgnvZ1tUaATfI/bvZOrQBgfIQzAACAARDOAAAABkA4A3bk1DgAgMMjnAEADJA/kMH8CGcAwOAJKsAcCGcAAAADIJwBAHfiXSqAzRDOAAAABkA4gwnx124AgPESzgAAAAZAOAMAABgA4QzYlVMlAdg0xyLmQjgD2CAvOACAE4QzAACAARDOOHTeKQAAgLsSzgAAAAZAOGPjvJMGAADCGQCwZgf9I5w/3gFzI5wxKw709ME4AmDoHKvGSTgDPsdEzl4YJ8BhMucwJ8IZQMcLAPpiLAHcmXlxb4QzAGBHXlABHB7hDNiWF2QADI1jE1MnnAHEAR8A2DzhDIDPEVLvoBZwZ3vZJ+w3sBrhDGBivDgCgHESzoBJE1QA1meTc6z5nSkSzgCAUfBiHJg64QwAAGAAhDMAAA7MO5qbZxtMh3DG7JjAAO7K3AiwecIZcFJetAEArJ9wBgAAMADCGTBJ3u0DAMZGOIMJEESgP/YnODz2t/VS3/ERzmBFJj44mMPed+yrMF/2f8ZCOANmw8EZABgy4QyYFQENYL02Mc+a25kK4QxgH7wA2D81g/Hauv/udX+23w+T7TJ8whmTZPIB+mAugdXYh2B/hDPYMAcuhsrYBKAPjid7J5zN3KZ3lk2vH4D1MccD7I9wxl04mAIAq/J6AvbvpOGsql5WVbdW1ZVLy55XVTdV1RXd1xPX20x2M8fJb459BqbDHNavIdVzSG2BOZjaPreXd85enuSCbZa/sLV2fvf1pn6bBQAAzNXUQtdenTSctdbeleQTh9AWmIWhTjZDbddWe2nnieccpE9jqQMAMD2rfObsWVX1/u60x9N3elJVXVxVl1XVZbfddtsKq4PdeVEN42F/BZiuneb4Veb+uRw3DhrOfj7JQ5Kcn+TmJC/Y6YmttUtba8dba8ePHDlywNXBHeaycwIAzM3cX+cdKJy11m5prX22tXZ7kpckeWS/zQLGZGwT6djau1er9muqdQEYCvPsauZQvwOFs6o6unT3m5NcudNz2Yw5DF4YAvsaJxgLAIdjyvPtXi6l/+okv5fkoVV1Y1U9I8lPVNUHqur9Sb4uyfeuuZ1M0JR3LGBc5j4fraP/c6/pugy5rn237aC/b8g12gT1GJe9XK3xwtba0dba3VtrZ7fWXtpae1pr7Staaw9rrX1Ta+3mw2gsMH3rOog4ONEH42iebPdxsJ36o5abs8rVGoGRMMkO3/I2sr2YKmMbYHfCGXAX+3kB5cXWOKxzOxkDdxh7LcbeftgL4/xg1O1wCGcAAAyGEMCcx4BwBhyaOU+2h2Vrjfusue03L7Y3DI/9cvqEMybLBAabZR8EDmrM88eY287mCWewRybb8ZrKZcL7WqexDADDJJwdgim/EJpy3zg52585Md5hmOybTIlwNnJzmZDm0k+YK/v46tRwfmxzmB7hbMleJ7mTPW+Tk6WJetz6GlvGAcCwbGpedjw4HOpMX4QzgJnzomJ36nO41BuYM+GMUXHQhnk4zH3dvAIwbHOap4WzGZjTgF4H9YO7sl8AU3BYc9mm58xNr5+9E84mbko745T6spu59HPuttvOfW/7sYylsbSTYTgxXsYwbsbQxt2Mvf0wRsIZjNxUDp5T6QebM4UxNIU+nMwc+pgcrJ9zqc2UTW0bTq0/YyCcrcigPRh1ow+HPY6MW2AMhj5XDb19jNvYx5dwtoOxb9h12W9dplbHqfUH9srYHybbhakzxvs3pJqu0pYh9aNPwtkh6WMAjXUQjrXdAGNn/p0325+pm+IYF84GagqDbch9GPo/Ax3TB95Zv2OXvHEwY2Es/wh90+uHdTCuYfqEM2DUvFgZnylus936NKRwzd7YXkzdXsa4/WAzhLMJ8q4L62RcMUbGLcyDC0UxdsLZgM3lHyNuSp/9nmsN+zbkOs79YjiQGNdz7v86+z7nusJWwhmTYGJnjuY07tfV1znVcK/mVJM59ZX1M54Ox9TrLJwBd7KOSW/qEykAsD2vAfZHOANYEwckAGA/hLND5IXaeoy5rie7wtuUDLU/c9oGjOdfAbB3ttUw2A7bG9MFSmzDYRDOAIBZ86J02qa2fafWn2SafToo4Qw4EBPpfLhyLKuwXcfPNmQ3xke/hDNgdBwIYP/8M2yA4RPOtnDg2owx1X1IbR1SWzg423F3e6nPWGs41navaq79hsS/BmF3wtkK7AQwTXPZt6fQzyH0YQhtABg7c+mCcAYxIcAUDXW/Hmq76I9tPA5z204+PzwOwhnAhhz0AObAx/IYMB7GxzZjO8YFiXDGkv1OCjs9v8/JxUTFVBjLMExb90376slNoUaH2Ycp1IvDI5zBIZnD5CyYw/TH7lj7N9Z2MwzGD4dFOAOYOC8qgGXmBA6DcXYwwhnswRAnmCG26aBW7cumajGlbQAMU18fOWB1asthEM5mwoQyHyfb1lMdC1PtF8zFHPfhOfb5sEy9tj5GMF3C2QGNcSCPsc2sznZfv03U2Ha9w2HUQr2Bw2K+mTfhbBvrPIXADsdBGTvrs11t1Xs91l1X242DGtvYGVt7YSfrGMtj3j+EM5LM95z2qfRjTNScKRnaeD7s9gyt/5ugBkyNMb1ZwhkcgqH+H52htIM72CaHYwh13ksbhtDOTZhrvxkX4/Rg1G13whl7YkfaPx/WXd1c+30Y+q7tkLbVkNqyTnPp56rmWqe59hvGTjhjUhyMNkftYd5cGGcc1AyG7aThrKpeVlW3VtWVS8vOqKq3VtU13ffT19vMeTBhwsLyvmC/2J367N/cajbW/h603evq72FdzGan9Yx1OwL7s5d3zl6e5IItyy5J8vbW2rlJ3t7dZ59MtPTJeAIYD3M2e2GczM9Jw1lr7V1JPrFl8ZOSvKK7/YokT+63WczVdu+YrDIxuWw3TI/9DmBz/Aup9TroZ87ObK3d3N3+WJIze2rPoAw5FEyRmjEVUxrLU+oLw2SMwXjYX9dv5QuCtNZakrbT41V1cVVdVlWX3XbbbauubnT6GsRD2hmG1BaATZnzXHiQvh9mvTa1bXxejN0YB+zFQcPZLVV1NEm677fu9MTW2qWtteOtteNHjhw54OoAGAP/QoI5GOLYHGKbDsOU/y0I83TQcPb6JBd1ty9K8pv9NIepmOPkNsc+A/vjH09zEKuOCWNqmmzXadrLpfRfneT3kjy0qm6sqmckeX6Sr6+qa5L8g+7+rNlBpsl2Zer6uPAOTPEU/sMyxz6zN8bGPJ16sie01i7c4aHH9dwWAIBeHLvkjbn++f9w7evYBC/ap8c25YSVLwjCfJg4AObFvD9ftj1shnA2EJuYBE28jNkYx+8Y2wyMx37nmGOXvNG8BAMjnMGGbPKA6GAM9oMhOaxt0cd6jJtxG9L2G1JbGA7hjMkx2cFm2Qf7oY4A8yOcDZADMkNnjO6PerEK4+euDqMmc6v73Pq7HyerjdqpQZ+Es5Mw2IBVzGkOWXdf51RL+rWOsWM8AusgnHXmNsnOrb8Amza1eXdq/dmkudRyt37OpQYnzK2/7J1wBgAAnTkHp7H3feztT4SzUZrCwGPzXC2SMTJ2DkbdwH7AOAhnMBMOStNgOwJ9MJfAMAlnADAzfb4w9yIfoD/C2T44AAGwKX0dgxzLYFyGuM8OsU1TIZwBAHDovMCHuxLOmD0HBwBgmdcGbIpwRq92mszGMMmNoY0A2zF/Qf/sV2yCcDZiJg3GYMyBHYB5O+xj1ZiOjWNq65gIZ/u0joFocDNVxjYnY4wAY2TuYl2Es5ExGQyb7QP9mtM+Nae+ArA94WyNHGgZGmMSYMF8CHewPwyHcAbbMEkBY2X+gnmZ0j4/pb4clHAGAAAjJMxMj3A2MHYyAOjHfo+pjsHApgln7MqBCgCmzbEehkM4GykTKQDQJ68t2C9jpn/CGb2Z2g46tf4A42U+ApgH4QzWwAup+bLtAWDzxno8Fs44qbEObgAAGBPhDACAQ+EPvnujTvMlnAHADHnxBzA8whkMnBdQAADzIJxlsy9+vfAGAAAS4QzomT84AAAcjHAGAMCh8oc82J5wBgAAMADCGfTMXwMBADgI4QwAAGAAhDMAAIABEM4AAAAGQDgDgIHx2VWAeRLOAAAABkA4AwAAGADhDAAAYACEMwAAgAEQzgAAAAbg1FV+uKquT/LpJJ9N8pnW2vE+GgUAADA3K4Wzzte11j7ew+8BAACYLac1cif+tw4AAGzGquGsJXlLVV1eVRf30SAAAIA5WvW0xq9urd1UVV+U5K1V9UettXctP6ELbRcnyQMe8IAVVwcAADBNK71z1lq7qft+a5LXJXnkNs+5tLV2vLV2/MiRI6usDgAAYLIOHM6q6vOr6rQTt5M8PsmVfTUMAABgTlY5rfHMJK+rqhO/51Wttd/upVUAAAAzc+Bw1lr7SJK/02NbAAAAZsul9AEAAAZAOAMAABgA4QwAAGAAhDMAAIABEM4AAAAGQDgDAAAYAOEMAABgAIQzAACAARDOAAAABkA4AwAAGADhDAAAYACEMwAAgAEQzgAAAAZAOAMAABgA4QwAAGAAhDMAAIABEM4AAAAGQDgDAAAYAOEMAABgAIQzAACAARDOAAAABkA4AwAAGADhDAAAYACEMwAAgAEQzgAAAAZAOAMAABgA4QwAAGAAhDMAAIABEM4AAAAGQDgDAAAYAOEMAABgAIQzAACAARDOAAAABkA4AwAAGADhDAAAYACEMwAAgAEQzgAAAAZAOAMAABgA4QwAAGAAhDMAAIABEM4AAAAGQDgDAAAYAOEMAABgAIQzAACAAVgpnFXVBVV1dVVdW1WX9NUoAACAuTlwOKuqU5L8XJInJDkvyYVVdV5fDQMAAJiTVd45e2SSa1trH2mt/XWS1yR5Uj/NAgAAmJdVwtlZSW5Yun9jtwwAAIB9qtbawX6w6ilJLmitPbO7/7Qkj2qtPWvL8y5OcnF396FJrj54c9fmfkk+vulGTIA69kMd+6GO/VDHfqhjP9SxH+rYD3Xsxxzr+MDW2pHtHjh1hV96U5Jzlu6f3S27k9bapUkuXWE9a1dVl7XWjm+6HWOnjv1Qx36oYz/UsR/q2A917Ic69kMd+6GOd7bKaY1/kOTcqnpQVd0jyVOTvL6fZgEAAMzLgd85a619pqqeleTNSU5J8rLW2gd7axkAAMCMrHJaY1prb0rypp7askmDPu1yRNSxH+rYD3Xshzr2Qx37oY79UMd+qGM/1HHJgS8IAgAAQH9W+cwZAAAAPZl9OKuqC6rq6qq6tqou2XR7NqWqrq+qD1TVFVV1WbfsjKp6a1Vd030/vVteVfWfupq9v6oesfR7Luqef01VXbS0/O92v//a7mdrt3WMRVW9rKpuraorl5ZtrG67rWPIdqjj86rqpm5MXlFVT1x67LldH6+uqm9YWr7t/txduOjd3fJf6S5ilKq6Z3f/2u7xYydbx5BV1TlV9c6q+lBVfbCqvqdbbkzuwy51NCb3oaruVVXvqar3dXX8t93y3vreZ32Hapc6vryqrlsaj+d3y+3Xu6iqU6rqvVX1hu6+8XgA29TReOxLa222X1lcyOTDSR6c5B5J3pfkvE23a0O1uD7J/bYs+4kkl3S3L0ny493tJyb5rSSV5NFJ3t0tPyPJR7rvp3e3T+8ee0/33Op+9gm7rWMsX0m+Jskjklw5hLrttI6hf+1Qx+cl+b5tnntet6/eM8mDun34lN325yS/muSp3e0XJfnO7vZ3JXlRd/upSX5lt3Vsuk57qOPRJI/obp+W5I+7vhiT/dTRmNxfHSvJfbrbd0/y7m4M9NL3Pus75K9d6vjyJE/Z5vn2693r+a+TvCrJGw4yVuY+Hnepo/HYV2033YAND6yvSvLmpfvPTfLcTbdrQ7W4PncNZ1cnOdrdPprk6u72i5NcuPV5SS5M8uKl5S/ulh1N8kdLyz/3vJ3WMaavJMdy51CxsbrttI5N1+iAdXxetn8hfKf9NIsrxn7VTvtzN1l/PMmp3fLPPe/Ez3a3T+2eVzutY9M1OkBNfzPJ1xuTvdXRmDx4DT8vyR8meVRffe+zvpuuzwHr+PJs/2LYfr1z/c5O8vYkj03yhj73xTmNx6117JYZjz19zf20xrOS3LB0/8Zu2Ry1JG+pqsur6uJu2ZmttZu72x9LcmZ3e6e67bb8xm2W77aOMdtk3aY2pp/VnaLwsrrjlNf91vFvJflka+0zW5bf6Xd1j/959/zR17E7PebhWfyV3Zg8oC11TIzJfelOfboiya1J3prFOwt99b3P+g7a1jq21k6Mxx/txuMLq+qe3TL79c5+Osn3J7m9u9/nvjib8Zi71vEE47EHcw9n3OGrW2uPSPKEJN9dVV+z/GBb/CmirbMBh7GOw6ZuK/n5JA9Jcn6Sm5O8YKOtGZGquk+S1yZ5dmvtU8uPGZN7t00djcl9aq19trV2fhZ/aX9kki/ZbIvGaWsdq+rLs3hX5kuSfGUWp4b9wJrbMOr9uqq+McmtrbXLN92WMduljsZjT+Yezm5Kcs7S/bO7ZbPTWrup+35rktdlcRC9paqOJkn3/dbu6TvVbbflZ2+zPLusY8w2WbfJjOnW2i3dC5Lbk7wkizGZ7L+O/zvJfavq1C3L7/S7use/sHv+aOtYVXfPIlC8srX2G91iY3KftqujMXlwrbVPJnlnFqd09dX3Pus7Ckt1vKC1dnNb+Kskv5CDj8e57NePSfJNVXV9ktdkcUrez8R43K+71LGqftl47M/cw9kfJDm3u4rOPbL4MObrN9ymQ1dVn19Vp524neTxSa7MohYXdU+7KIvPXaRb/u3d1XEeneTPu7eZ35zk8VV1ene6z+OzOK/65iSfqqpHd1fc+fYtv2u7dYzZJuu20zpG58QE3PnmLMZksujjU2txlasHJTk3iw8Pb7s/d39de2eSp3Q/v7VeJ+r4lCTv6J6/0zoGrRsnL01yVWvtp5YeMib3Yac6GpP7U1VHquq+3e17Z/G5vavSX9/7rO9g7VDHP1p6kVpJnpw7j0f79Rattee21s5urR3LYqy8o7X2rTEe92WHOn6b8dijg35YbSpfWVzh5Y+zOA/+hzbdng3V4MFZXFXofUk+eKIOWZz3/PYk1yR5W5IzuuWV5Oe6mn0gyfGl3/UdSa7tvp6+tPx4Fjvqh5P8bO74B+jbrmMsX0lencXpTX+TxTnOz9hk3XZbx5C/dqjjL3V9eH8WE+/Rpef/UNfHq9Ndxalbvu3+3I3x93T1/bUk9+yW36u7f233+INPto4hfyX56ixO83h/kiu6rycak73V0ZjcXx0fluS9Xb2uTPLDffe9z/oO9WuXOr6jG49XJvnl3HFFR/v1yWv6tbnjQhbGYz91NB57+jrRWQAAADZo7qc1AgAADIJwBgAAMADCGQAAwAAIZwAAAAMgnAEAAAyAcAYAADAAwhkAAMAACGcAAAAD8P8BejRKKkFFCi0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1080x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean_cov=12.01, breadth_cov=0.9461779111683469\n",
      "Indel stats\n",
      "2725859 indels of length 1\n",
      "272869 indels of length 2\n",
      "47507 indels of length 3\n",
      "14045 indels of length 4\n",
      "5897 indels of length 5\n",
      "2866 indels of length 6\n",
      "1615 indels of length 7\n",
      "1016 indels of length 8\n",
      "650 indels of length 9\n",
      "424 indels of length 10\n",
      "303 indels of length 11\n",
      "227 indels of length 12\n",
      "191 indels of length 13\n",
      "143 indels of length 14\n",
      "86 indels of length 15\n",
      "61 indels of length 16\n",
      "35 indels of length 17\n",
      "48 indels of length 18\n",
      "27 indels of length 19\n",
      "26 indels of length 20\n",
      "22 indels of length 21\n",
      "18 indels of length 22\n",
      "17 indels of length 23\n",
      "14 indels of length 24\n",
      "22 indels of length 25\n",
      "10 indels of length 26\n",
      "13 indels of length 27\n",
      "17 indels of length 28\n",
      "6 indels of length 29\n",
      "6 indels of length 30\n",
      "9 indels of length 31\n",
      "2 indels of length 32\n",
      "6 indels of length 33\n",
      "7 indels of length 34\n",
      "5 indels of length 35\n",
      "5 indels of length 36\n",
      "3 indels of length 37\n",
      "2 indels of length 38\n",
      "3 indels of length 39\n",
      "7 indels of length 40\n",
      "2 indels of length 41\n",
      "2 indels of length 42\n",
      "2 indels of length 43\n",
      "1 indels of length 44\n",
      "1 indels of length 45\n",
      "3 indels of length 46\n",
      "2 indels of length 47\n",
      "2 indels of length 48\n",
      "5 indels of length 49\n",
      "3 indels of length 51\n",
      "2 indels of length 52\n",
      "4 indels of length 53\n",
      "2 indels of length 59\n",
      "1 indels of length 60\n",
      "2 indels of length 62\n",
      "2 indels of length 68\n",
      "1 indels of length 69\n",
      "1 indels of length 70\n",
      "2 indels of length 75\n",
      "1 indels of length 80\n",
      "1 indels of length 81\n",
      "2 indels of length 85\n",
      "1 indels of length 87\n",
      "1 indels of length 90\n",
      "1 indels of length 92\n",
      "1 indels of length 98\n",
      "1 indels of length 99\n",
      "1 indels of length 106\n",
      "1 indels of length 109\n",
      "1 indels of length 115\n",
      "1 indels of length 121\n",
      "1 indels of length 141\n",
      "1 indels of length 145\n",
      "1 indels of length 146\n",
      "1 indels of length 157\n",
      "1 indels of length 175\n",
      "1 indels of length 187\n",
      "1 indels of length 206\n",
      "1 indels of length 283\n",
      "Indels near homopolymers: 33%\n"
     ]
    }
   ],
   "source": [
    "print(\"========== pacbio_40x.fq.gz  ==========\")\n",
    "\n",
    "\n",
    "samfile = pysam.AlignmentFile(\"pacbio_sorted.bam\", \"rb\")\n",
    "genome = next(SeqIO.parse(\"reference.fasta\", \"fasta\")).seq\n",
    "\n",
    "\n",
    "# coverage & errors\n",
    "cov = np.zeros(shape=(len(genome)))\n",
    "SNP = pd.DataFrame(0, index=['A', 'C', 'G', 'T', '-'], columns=['A', 'C', 'G', 'T', '-'])\n",
    "quals = pd.DataFrame(0, index=['A', 'C', 'G', 'T', '-'], columns=['A', 'C', 'G', 'T', '-'])\n",
    "\n",
    "#print(cov)\n",
    "#print(SNP)\n",
    "\n",
    "delta = 0\n",
    "indels = Counter() # length -> count\n",
    "homopoly = 0\n",
    "\n",
    "for i, x in enumerate(samfile.pileup('gi|48994873|gb|U00096.2|', fastafile=pysam.FastaFile('reference.fasta'))):\n",
    "    pos = x.reference_pos\n",
    "    cov[i] = x.get_num_aligned()\n",
    "    gen = genome[pos].upper()\n",
    "    for base in x.get_query_sequences(add_indels=True):\n",
    "        j = base[0].upper()\n",
    "        if j != gen and j in \"AGCT-\":\n",
    "            SNP.loc[j, gen] += 1\n",
    "        INS = True\n",
    "        if len(base) > 1:\n",
    "            if base[1] == '-':\n",
    "                INS = False\n",
    "            num = \"\"\n",
    "            for ind in base[2:]:\n",
    "                if ind in \"0123456789\":\n",
    "                    num += ind\n",
    "                else:\n",
    "                    break\n",
    "            num = int(num)\n",
    "            indels[num] += 1\n",
    "            \n",
    "            if gen == genome[pos-1].upper() or gen == genome[pos+1].upper() or \\\n",
    "            len(base[len(base)-num:]) > 1 and base[len(base)-num:][0] == base[len(base)-num:][1]:\n",
    "                homopoly += num\n",
    "            for ind in base[len(base)-num:]:\n",
    "                if INS:\n",
    "                    SNP.loc['-', ind.upper()] += 1\n",
    "                else:\n",
    "                    SNP.loc[ind.upper(), '-'] += 1\n",
    "\n",
    "print(\"Number of substitutions:\")\n",
    "print(SNP)\n",
    "total = sum(cov)\n",
    "print(f'mean_read_error={round(SNP.iloc[:4, :].values.sum()/total*100, 5)}%')\n",
    "\n",
    "\n",
    "covs = [np.mean(cov[i:i+1000]) for i in range(0, len(cov), 1000)]\n",
    "plt.figure(figsize=(15, 5))\n",
    "plt.bar(range(len(covs)), height=covs)\n",
    "plt.xticks(range(0, len(covs), 500), range(0, len(cov), 500000))\n",
    "plt.title(\"Mean coverage\")\n",
    "plt.show()\n",
    "print(f'mean_cov={round(np.mean(cov), 2)}, breadth_cov={np.count_nonzero(cov)/len(cov)}')\n",
    "\n",
    "print('Indel stats')\n",
    "for k, v in sorted(indels.items()):\n",
    "    print(f'{v} indels of length {k}')\n",
    "    \n",
    "print(f'Indels near homopolymers: {round(homopoly/sum(indels.values())*100)}%')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "silver-character",
   "metadata": {},
   "source": [
    "# Oxford Nanopores"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "loving-trading",
   "metadata": {},
   "source": [
    "##### >minimap2 -x map-ont -a -t 6 -o ONT_2D.sam reference.fasta ONT_2D.fastq.gz\n",
    "\n",
    "##### >samtools flagstat -@ 6 ONT_2D.sam > ONT_2D.stat\n",
    "\n",
    "32103 + 0 in total (QC-passed reads + QC-failed reads)\n",
    "\n",
    "1072 + 0 secondary\n",
    "\n",
    "221 + 0 supplementary\n",
    "\n",
    "0 + 0 duplicates\n",
    "\n",
    "29565 + 0 mapped (92.09% : N/A)\n",
    "\n",
    "0 + 0 paired in sequencing\n",
    "\n",
    "0 + 0 read1\n",
    "\n",
    "0 + 0 read2\n",
    "\n",
    "0 + 0 properly paired (N/A : N/A)\n",
    "\n",
    "0 + 0 with itself and mate mapped\n",
    "\n",
    "0 + 0 singletons (N/A : N/A)\n",
    "\n",
    "0 + 0 with mate mapped to a different chr\n",
    "\n",
    "0 + 0 with mate mapped to a different chr (mapQ>=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "respiratory-scott",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "========== ONT_2D.fastq.gz  ==========\n",
      "Number of substitutions:\n",
      "        A       C       G       T       -\n",
      "A       0    6938   29590    5831  648838\n",
      "C    5808       0    4271   46647  561305\n",
      "G   46844    4346       0    5623  560507\n",
      "T    5727   29410    6996       0  610449\n",
      "-  123633  112716  120460  111910       0\n",
      "mean_read_error=2.6327%\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2cAAAE/CAYAAADCCbvWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAhXElEQVR4nO3df7B0d10f8PeHJPxQUBJ5ZB6S4AMYwagY6GOAwVoNgoFxBFs6BRVThGaqMBVFJdQZxba26hSjrVYIBUHlhwoyUEQhQCy1VeCJJCEhYAIJkhjIgxCBalHCt3/seczmcn/v2bvnnH29Znbu7tnde77ncz7n7L7vnj23WmsBAABgte6y6gEAAAAgnAEAAAyCcAYAADAAwhkAAMAACGcAAAADIJwBAAAMgHAGAAAwAMIZAPtSVTdW1d9V1X02TH9vVbWqOrKioQHAKAlnACzihiRPPXGjqr4hyZesbjirU1Unr3oMAIybcAbAIn4zyffP3b4gyW/MP6Cq7lZV/7mq/qKqPl5VL6qqe3T3nVpVb6qq41X1qe76GXPP/aOq+vdV9b+r6jNV9daNn9RtmNcTq+qKqvp0VX2oqs7vpt+vqt5YVZ+squur6l/NTf/bqjpt7nc8rKo+UVWndLd/oKqu7cb3lqr6qrnHtqp6VlVdl+S6btovV9VHuzFcXlX/eO7x96iqV3S/69qq+omqumnu/vtV1eu6etxQVf9mb6sDgDETzgBYxJ8m+bKq+tqqOinJU5L81obH/FySr0lyTpKvTnJ6kp/q7rtLkl9P8lVJ7p/kb5P8yobnf0+Spyf5yiR3TfJjmw2kqs7NLBj+eJJ7J/mWJDd2d78myU1J7pfkyUn+Y1Wd11r7yyR/kuSfbZjfa1trf19VT0zyb5P80ySHkvyvJK/eMOsnJXlEkrO72+/plvW0JK9K8rtVdffuvp9OciTJA5M8Nsn3zY3/Lkn+R5Iruxo9Jslzquo7NlteAKanWmurHgMAI1RVNyZ5ZpJHJvnSJP8zyXOTPD7J3yd5QJKPJPlskoe21j7UPe9RSV7VWnvAJr/znCSXtdZO7W7/UZK3tdb+Q3f7h5J8V2vt/E2e++Ikf9Na+5EN08/MLKTdu7X2mW7af0pyuLX2L6vqmUm+p7V2XlVVkr9I8r2ttXdW1R9kFtRe2j3vLt3yfG1r7SNV1ZI8prX2jm3q9Kkk39pau7KqPpzkB1trb+nue2aSF7TWzqiqRyT53dba/eee+/wkX9Nae/pWvx+A6XB8PACL+s0k78wsjP3GhvsOZfYdtMtnuSdJUklOSpKq+pIkFyc5P8mp3f33qqqTWmu3d7c/Nvf7/ibJPbcYx5lJ3rzJ9Psl+eSJYNb5SJKj3fXXJfmvVXU4s0/4vpDZJ2TJ7BO9X66qF849tzL7ZOsj3e2Pzs+sqn4syTO6+bYkX5bkxKGY99vw+PnrX5XkflV129y0k+bGAsDECWcALKT7BOmGJE/ILJTM+0Rmhyp+XWvt5k2e/twkD07yiNbax7pPzt6bWQDaq48medAm0/8yyWlVda+5gHb/JDd34/9UVb01yb9I8rVJXtPuOKzko0l+trX2ym3m+w+HoHTfL/uJzA5JvKa19oXuk7MTy3NLkjOSvL+7feaG8d/QWjtrV0sLwOT4zhkAfXhGkvNaa/93fmJr7QtJXpLk4qr6yiSpqtPnvkd1r8zC223dSTl+eoExvDTJ06vqMVV1l24+D2mtfTTJ/0nyn6rq7lX10G6889+Ne1VmJzZ5cnf9hBcleX5VfV039i+vqn++zRjuleTzSY4nObmqfiqzT85O+J3u951aVacnefbcfe9O8pmqel534pCTqurrq+qb9lELAEZIOANgYa21D7XWjm1x9/OSXJ/kT6vq00neltmnZUnyS0nukdknbH+a5A8XGMO7MztxyMVJ/jqz78CdOLPiUzM7EcdfJnl9kp9urb1t7ulvTHJWko+11q6c+52vT/LzSV7Tjf3qzL5Tt5W3dMvw55kd9vj/cudDF/9dZicmuSGzOrw2yee6ed2e5DszO5nIDZnV5L8n+fJdFwGAUXNCEABYkar6wSRPaa39k1WPBYDV88kZAByQqjpcVY/uDrt8cGbfuXv9qscFwDA4IQgAHJy7JnlxZme2vC2z/7/231Y5IACGw2GNAAAAA+CwRgAAgAEQzgAAAAbgQL9zdp/73KcdOXLkIGcJAAAwGJdffvknWmuHNrvvQMPZkSNHcuzYVv8GBwAAYNqq6iNb3eewRgAAgAEQzgAAAAZAOAMAABgA4QwAAGAAhDMAAIABEM4AAAAGQDgDAAAYAOEMAABgAIQzAACAARDOAAAABkA4AwAAGADhDAAAYACEMwAAgAEQzgAAAAZAOAMAABgA4QwAAGAAhDMAAIABEM4AAAAGQDgDAAAYAOEMAABgAHYMZ1V196p6d1VdWVXXVNXPdNNfXlU3VNUV3eWcpY8WAABgok7exWM+l+S81tpnq+qUJH9cVX/Q3ffjrbXXLm94AAAA62HHcNZaa0k+2908pbu0ZQ4KAABg3ezqO2dVdVJVXZHk1iSXttbe1d31s1V1VVVdXFV3W9YgAQAApm5X4ay1dntr7ZwkZyQ5t6q+PsnzkzwkyTclOS3J8zZ7blVdWFXHqurY8ePH+xk1AADAxOzpbI2ttduSXJbk/NbaLW3mc0l+Pcm5Wzznktba0dba0UOHDi08YAAAgCnazdkaD1XVvbvr90jy2CQfqKrD3bRK8qQkVy9vmAAAANO2m7M1Hk7yiqo6KbMw9zuttTdV1Tuq6lCSSnJFkn+9vGECAABM227O1nhVkodtMv28pYwIAABgDe3pO2cAAAAsh3AGAAAwAMIZAADAAAhnADACRy76/VUPAYAlE84AAAAGQDgDAAAYAOEMAABgAIQzAACAARDOAAAABkA4AwAAGADhDAbKabMBANaLcAYAADAAwhkAAMAACGcAAAADIJwBAAAMgHAGwNpy4h0AhkQ4AwAAGADhDAAAYACEMwAAgAEQzibAdyYAAGD8hDMAAIABEM4AAAAGQDgDAAAYAOEMAABgAISzjpNqAAAAqyScAQAADMCO4ayq7l5V766qK6vqmqr6mW76A6rqXVV1fVX9dlXddfnDBQAAmKbdfHL2uSTntda+Mck5Sc6vqkcm+fkkF7fWvjrJp5I8Y2mjBAAAmLgdw1mb+Wx385Tu0pKcl+S13fRXJHnSMgYIAACwDnb1nbOqOqmqrkhya5JLk3woyW2ttc93D7kpyelLGSEAAMAa2FU4a63d3lo7J8kZSc5N8pDdzqCqLqyqY1V17Pjx4/sbJQAAwMTt6WyNrbXbklyW5FFJ7l1VJ3d3nZHk5i2ec0lr7Whr7eihQ4cWGSsAE+XfmQDA7s7WeKiq7t1dv0eSxya5NrOQ9uTuYRckecOSxggr4c0iMBT2RwDrYTefnB1OcllVXZXkPUkuba29KcnzkvxoVV2f5CuSvHR5w+SgeAMA68G2DoyJfRbr4uSdHtBauyrJwzaZ/uHMvn8GAADAgvb0nTMAAACWQzgDAIARcHjn9AlnAAAAAyCcAQAADIBwBgAwYA5lYzv6Y1qEMxi5Me2UxzRWgHVh3wzDIZzBmvEizEHa2G/6DwC2JpwxSd4AMmX6G5bPdgasgnAGe7SqF2xvFKbHOoWDYVsDxkI4A1bCmyWGRD8CMATCGQDABPmjA4yPcAYAADAAwhkAk+BTAoDls69dLuGMhdhAAQCgH8LZJgQOWE/rsu2vy3ICwNgIZ7AC3hzDatkGARgi4YzeedMDjIF9FTA19mvjJ5wBrMDUXkCHsjxDGQcsm17fnvpM0zqsV+GML7IOjb9basEY6Vv6sC59tC7Lyeas/62NpTZjGeduCWeQ6W3YAACMj3AGPVjncLfMZV/0d6/zemE6FunjvrcB2xTAcglnwFrwphIOlm0OYO+EM3ZtWS+0XsABgGXxPoMxEc5gAsbwwjOGMTJ9+nD9WOf9Uk9YLuEM1pAXV2Co7J/uTD04CGPoszGMsQ/CGexgXXYGjIu+5KD02Wv6lnWgz1nEjuGsqs6sqsuq6v1VdU1V/XA3/QVVdXNVXdFdnrD84cLBsGMFAOCg7eaTs88neW5r7ewkj0zyrKo6u7vv4tbaOd3lzUsbJXciOAzTGNbLGMa4nbGPH4Bp8zrFonYMZ621W1prf9Zd/0ySa5OcvuyBAQAwPAIILM+evnNWVUeSPCzJu7pJz66qq6rqZVV1at+DY73Y2bMK+g5g+uzr18MU1vOuw1lV3TPJ65I8p7X26SS/luRBSc5JckuSF27xvAur6lhVHTt+/PjiI2aSprAxsRx6YzyGvq6GPj5gnIaybxnKOFjMrsJZVZ2SWTB7ZWvt95Kktfbx1trtrbUvJHlJknM3e25r7ZLW2tHW2tFDhw71NW7onZ0aq6DvmDL9DbA3uzlbYyV5aZJrW2u/ODf98NzDvjvJ1f0PDwDGQxhhlfbTf3p2fIa0zvYyliGNe8h288nZo5M8Lcl5G06b/wtV9b6quirJtyX5kWUOdIw0IfP0AzB0q95PrXr+bM26Yd58PxxEb6xT/+3mbI1/3Fqr1tpD50+b31p7WmvtG7rp39Vau+UgBgx7sU4bM7A/9hPjd9BvFAGWZU9na6R/XkRg2tZ5G1/nZWdr+gJYprHvY4QzYFtj38kNyVa1XGaNx7T+xjTWzYx9/PRHLwD7JZyNmJ0/Q6EXWTU9CONhe909tVo/whnAHnmxJFl+H+iz4XOmummyrlgl4awHNuIZdWARq+yfMfTuGMYIcMKQ91lDHhsIZ+zK0Hdkm41v6GMGtrbu2+8iy7/utQOGzT5qe8LZGrNxjId1NTzWCay3Ze0D7FtgvQlnG9gpAuvE4aS7t914x7YswHo7yH3WMuY15X2ucAYTM+UdFgytv4c2nlU56DoMqe5DGgswM+btUjgbKH9loE/WPRwc29vO1Iix0Ku703ed1rnuwhmDtY4b5jouM6un7wBm7A9ZNeFsn2y87FdfvaMHOUiLfN9Krx4ctWZspt6zU1++3VKH3RPOoGdj/pLtQYx9ijvoKS7T2FknrDP9vzlfGRmHda+pcDYS696o66qP9a53gLGy/xq3Vf2/Pn2zP+o2DMLZyNmQDsZe6+z/3wA7sT0DsJFwdoAWfSH2Qs4y6S+GRk9uT33GzfqbljGvz/2OfczLPGTC2Q42azzNuDMnvQCYLvtmgOUQzoBBGdubvrGMdz/jPKhlG0sNgcXsdlsfylcJ9moo42DchLOJsoNgnSzrBR/WyTpsH0Nfxj7HN/RlHQM1ZBWEM9aaHe/uqNO4WX8Aw7HKfbJ/mTN8wtk2xvz/qpgOvQH7s9W2M9ZDphiGVfXDUPtwaOMa2nhgr4Qz7MgGZkjrY0hjWcReluPIRb8/meXmDgIZTNsyt1ln2+YgCWd74B8CMyV6kd0awiE4qzqSwXayHjau5ymt9ymc2GdK62PeVJeLxQhne2RDgmFwEpDh2k3Np7ReprQsG0152YZGrQ/GOtd5qMs+1HGtinDGjhbZaGxwd6Ye7Mey+kY/Anux0z5j6vuUqS8fwyCcsW9T2UkNcTmGOCZYZ7bJ8TuoYLGfw2L1F1O11+98I5zBytgJwXpuB+u4zMswxE+Ulx3GpnoW6XV4Az+kcQ9pLHyxHcNZVZ1ZVZdV1fur6pqq+uFu+mlVdWlVXdf9PHX5w6VvYziD2RBfgNed2gHA/nkd3Zy67O6Ts88neW5r7ewkj0zyrKo6O8lFSd7eWjsrydu72/Ro7A069vEDwDJN8XVyiss0RMs+O6b1uDo7hrPW2i2ttT/rrn8mybVJTk/yxCSv6B72iiRPWtIY2SMbFKug78bN+oPlsX1Bf6a+Pe3pO2dVdSTJw5K8K8l9W2u3dHd9LMl9t3jOhVV1rKqOHT9+fJGxsker+P9ATMOQT4Wun1drneq/Tss6NGoPrKtdh7OqumeS1yV5Tmvt0/P3tdZakrbZ81prl7TWjrbWjh46dGihwcKivOD3Sz2Hxz+CZSdjXo9jHjv90w9M0a7CWVWdklkwe2Vr7fe6yR+vqsPd/YeT3LqcIU6DHcjBUWsSfQAn2BaYN6R+GNJYYCh2c7bGSvLSJNe21n5x7q43Jrmgu35Bkjf0P7zV8Y+Xx2cIdd9qDEMYGwfLOofpmcp2PZXl2Giqy8Xmprq+d/PJ2aOTPC3JeVV1RXd5QpKfS/LYqrouybd3t9mDg2iqqTYu621sfT3k7/CNyVhrtOzv/461LnvV93L6IywwRLs5W+Mft9aqtfbQ1to53eXNrbW/aq09prV2Vmvt21trnzyIATNtXvAYOj3KEOlL9mqzntFH02Odjs+eztYIrI8+d+g7/S5/wYa90fcA0yScAYyUfxQKB8N2tjcO4aUP67q+hTOYkO12ZOu6k4M+DXk7OqixLWM+Q64r62kMPbnMo1IOyhjGeNCEs544LGv51Imx8SaWoVjm94v6OuGJ7QVAOLuTve7EfWzPdqxHVm0VPajvARiCsb4eCWcj5S+MjJleA/oy5sM5h2TqyzcW1gPCGXt2kGfxm7p1X/6NVlmPqa2L/SzP1GoAe2UbAFZNOFtTB/0CNJUXvKksB/tj/QMwJF6Xpkc4AwbHiw2wiHXchzhbL2OgF3cmnDEKNmYAoE/eWzBEwtmSTH2Dn/ryMW76c71Z/wCMlXAGE+ENKQdJv43PqtfZqud/0NZteRelXjAjnB0AOxwYP9sxfdBHAGxHOJsYL/zsxbL75SD7cbfzso3sjXoBU2TfxlAJZ0s29o1/7ONfV9Ybq6YHAWDvhLMFeQOy3qa4/qe4TGMxhNoPYQzrZLN6WwcsQv/AuAlnABwobx5hd1Yd3m2r42FdTYdwBqydZb6IeYFkneh3YDP2DfsnnNErGyPAzuwrWYT+ucNUa9Hnck21Rsk0l004AzgAfbyATPFFaD/UAeDg2fceDOEM6I0dNwBsbpWvkV6fx0M4A2BPvMizSlPrv6ktz0ZTX75lOFEztVtPwhnskp0ksF/73X+s835nnZd9SKwHOFjCGSthZw/9s10BwLgJZyPjzdf0WKfDZd0AixryfqSvsQ15GWFshDNYkBelO1OPxakhuzWFXpnCMgD0ZcdwVlUvq6pbq+rquWkvqKqbq+qK7vKE5Q4TYD14o8pB2KnP9CHAauzmk7OXJzl/k+kXt9bO6S5v7ndYAMAqCGawPFttX7Y7TtgxnLXW3pnkkwcwFpZoKv9bw84LAICpWuQ7Z8+uqqu6wx5P3epBVXVhVR2rqmPHjx9fYHZMhYDFbugTTtiuF/RJP9SRg6LXxs36W779hrNfS/KgJOckuSXJC7d6YGvtktba0dba0UOHDu1zdgAA/Rr7G82xjx/4YvsKZ621j7fWbm+tfSHJS5Kc2++wDpad27hZfwDryf6fVdJ/LMO+wllVHZ67+d1Jrt7qsQDAevMmlo3G3hPLGP/Ya0I/dnMq/Vcn+ZMkD66qm6rqGUl+oareV1VXJfm2JD+y5HHCrtixDY91wjrS97A429G4WX/7c/JOD2itPXWTyS9dwlgAAADW1iJna6Rn/sIA42KbHQbrYTXUHaB/whkAAMAACGesJX/x3T+1g+U4ctHv276AQbFPOnjCGQAAbGNdQsq6LOeQCWcAAAADIJwxOov+VcdfhQAOjn0uwO4JZwAAAAMgnAEAAAyAcMbKONSFRS2zh/QnU6Sv2Q99AwdHOAMAABgA4QwAAGAAhDNGySEWAABMjXA2MEIHAACsJ+EMAABgAIQzAACAARDOAIDJ8TUBYIyEM4CeeVMI0A/7UzVYN8IZAADAAAhnJPFXGQAAWDXhDAAAWBkfEtxBOAMAABgA4QwAAGAAhDMAAIABEM4AAAAGQDgDAAAYAOEMAABgAHYMZ1X1sqq6taqunpt2WlVdWlXXdT9PXe4wAQAApm03n5y9PMn5G6ZdlOTtrbWzkry9uw0AAMA+7RjOWmvvTPLJDZOfmOQV3fVXJHlSv8MCAABYL/v9ztl9W2u3dNc/luS+PY0HAABgLS18QpDWWkvStrq/qi6sqmNVdez48eOLzg4AAGCS9hvOPl5Vh5Ok+3nrVg9srV3SWjvaWjt66NChfc4OAABg2vYbzt6Y5ILu+gVJ3tDPcAAAANbTbk6l/+okf5LkwVV1U1U9I8nPJXlsVV2X5Nu72wAAAOzTyTs9oLX21C3uekzPYwEAAFhbC58QBAAAgMUJZwAAAAMgnAEAAAyAcAYAADAAwhkAAMAACGcAAAADIJwBAAAMgHAGAAAwAMIZAADAAAhnAAAAAyCcAQAADIBwBgAAMADCGQAAwAAIZwAAAAMgnAEAAAyAcAYAADAAwhkAAMAACGcAAAADIJwBAAAMgHAGAAAwAMIZAADAAAhnAAAAAyCcAQAADIBwBgAAMADCGQAAwAAIZwAAAANw8iJPrqobk3wmye1JPt9aO9rHoAAAANbNQuGs822ttU/08HsAAADWlsMaAQAABmDRcNaSvLWqLq+qC/sYEAAAwDpa9LDGb26t3VxVX5nk0qr6QGvtnfMP6ELbhUly//vff8HZAQAATNNCn5y11m7uft6a5PVJzt3kMZe01o621o4eOnRokdkBAABM1r7DWVV9aVXd68T1JI9LcnVfAwMAAFgnixzWeN8kr6+qE7/nVa21P+xlVAAAAGtm3+GstfbhJN/Y41gAAADWllPpAwAADIBwBgAAMADCGQAAwAAIZwAAAAMgnAEAAAyAcAYAADAAwhkAAMAACGcAAAADIJwBAAAMgHAGAAAwAMIZAADAAAhnAAAAAyCcAQAADIBwBgAAMADCGQAAwAAIZwAAAAMgnAEAAAyAcAYAADAAwhkAAMAACGcAAAADIJwBAAAMgHAGAAAwAMIZAADAAAhnAAAAAyCcAQAADIBwBgAAMAALhbOqOr+qPlhV11fVRX0NCgAAYN3sO5xV1UlJfjXJ45OcneSpVXV2XwMDAABYJ4t8cnZukutbax9urf1dktckeWI/wwIAAFgvi4Sz05N8dO72Td00AAAA9qhaa/t7YtWTk5zfWntmd/tpSR7RWnv2hsddmOTC7uaDk3xw/8Ndmvsk+cSqBzEB6tgPdeyHOvZDHfuhjv1Qx36oYz/UsR/rWMevaq0d2uyOkxf4pTcnOXPu9hndtDtprV2S5JIF5rN0VXWstXZ01eMYO3Xshzr2Qx37oY79UMd+qGM/1LEf6tgPdbyzRQ5rfE+Ss6rqAVV11yRPSfLGfoYFAACwXvb9yVlr7fNV9ewkb0lyUpKXtdau6W1kAAAAa2SRwxrTWntzkjf3NJZVGvRhlyOijv1Qx36oYz/UsR/q2A917Ic69kMd+6GOc/Z9QhAAAAD6s8h3zgAAAOjJ2oezqjq/qj5YVddX1UWrHs+qVNWNVfW+qrqiqo51006rqkur6rru56nd9Kqq/9LV7Kqqevjc77mge/x1VXXB3PR/1P3+67vn1nbzGIuqellV3VpVV89NW1ndtpvHkG1RxxdU1c1dT15RVU+Yu+/53TJ+sKq+Y276pttzd+Kid3XTf7s7iVGq6m7d7eu7+4/sNI8hq6ozq+qyqnp/VV1TVT/cTdeTe7BNHfXkHlTV3avq3VV1ZVfHn+mm97bsfdZ3qLap48ur6oa5fjynm2673kZVnVRV762qN3W39eM+bFJH/diX1traXjI7kcmHkjwwyV2TXJnk7FWPa0W1uDHJfTZM+4UkF3XXL0ry8931JyT5gySV5JFJ3tVNPy3Jh7ufp3bXT+3ue3f32Oqe+/jt5jGWS5JvSfLwJFcPoW5bzWPoly3q+IIkP7bJY8/uttW7JXlAtw2ftN32nOR3kjylu/6iJD/YXf+hJC/qrj8lyW9vN49V12kXdTyc5OHd9Xsl+fNuWfRkP3XUk3urYyW5Z3f9lCTv6nqgl2Xvs75DvmxTx5cnefImj7ddb1/PH03yqiRv2k+vrHs/blNH/dhXbVc9gBU31qOSvGXu9vOTPH/V41pRLW7MF4ezDyY53F0/nOSD3fUXJ3nqxscleWqSF89Nf3E37XCSD8xN/4fHbTWPMV2SHMmdQ8XK6rbVPFZdo33W8QXZ/I3wnbbTzM4Y+6ittuduZ/2JJCd30//hcSee210/uXtcbTWPVddoHzV9Q5LH6sne6qgn91/DL0nyZ0ke0dey91nfVddnn3V8eTZ/M2y73rp+ZyR5e5Lzkrypz21xnfpxYx27afqxp8u6H9Z4epKPzt2+qZu2jlqSt1bV5VV1YTftvq21W7rrH0ty3+76VnXbbvpNm0zfbh5jtsq6Ta2nn90dovCyuuOQ173W8SuS3NZa+/yG6Xf6Xd39f909fvR17A6PeVhmf2XXk/u0oY6JntyT7tCnK5LcmuTSzD5Z6GvZ+6zvoG2sY2vtRD/+bNePF1fV3bpptuut/VKSn0jyhe52n9vi2vRjvriOJ+jHHqx7OOMO39xae3iSxyd5VlV9y/ydbfaniLbMARzEPA6aui3k15I8KMk5SW5J8sKVjmZEquqeSV6X5DmttU/P36cnd2+TOurJPWqt3d5aOyezv7Sfm+Qhqx3ROG2sY1V9fWafyjwkyTdldmjY85Y8hlFv11X1nUluba1dvuqxjNk2ddSPPVn3cHZzkjPnbp/RTVs7rbWbu5+3Jnl9Zi+iH6+qw0nS/by1e/hWddtu+hmbTM828xizVdZtMj3dWvt494bkC0lekllPJnuv418luXdVnbxh+p1+V3f/l3ePH20dq+qUzALFK1trv9dN1pN7tFkd9eT+tdZuS3JZZod09bXsfdZ3FObqeH5r7ZY287kkv5799+O6bNePTvJdVXVjktdkdkjeL0c/7tUX1bGqfks/9mfdw9l7kpzVnUXnrpl9GfONKx7TgauqL62qe524nuRxSa7OrBYXdA+7ILPvXaSb/v3d2XEemeSvu4+Z35LkcVV1ane4z+MyO676liSfrqpHdmfc+f4Nv2uzeYzZKuu21TxG58QOuPPdmfVkMlvGp9TsLFcPSHJWZl8e3nR77v66dlmSJ3fP31ivE3V8cpJ3dI/fah6D1vXJS5Nc21r7xbm79OQebFVHPbk3VXWoqu7dXb9HZt/buzb9LXuf9R2sLer4gbk3qZXkSblzP9quN2itPb+1dkZr7UhmvfKO1tr3Rj/uyRZ1/D792KP9flltKpfMzvDy55kdB/+Tqx7PimrwwMzOKnRlkmtO1CGz457fnuS6JG9Lclo3vZL8alez9yU5Ove7fiDJ9d3l6XPTj2a2oX4oya/kjn+Avuk8xnJJ8urMDm/6+8yOcX7GKuu23TyGfNmijr/ZLcNVme14D889/ie7ZfxgurM4ddM33Z67Hn93V9/fTXK3bvrdu9vXd/c/cKd5DPmS5JszO8zjqiRXdJcn6Mne6qgn91bHhyZ5b1evq5P8VN/L3md9h3rZpo7v6Prx6iS/lTvO6Gi73rmm35o7TmShH/upo37s6XJiYQEAAFihdT+sEQAAYBCEMwAAgAEQzgAAAAZAOAMAABgA4QwAAGAAhDMAAIABEM4AAAAGQDgDAAAYgP8Pll0k/WRYSaYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1080x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean_cov=21.11, breadth_cov=0.9968482706224036\n",
      "Indel stats\n",
      "1194620 indels of length 1\n",
      "381204 indels of length 2\n",
      "155145 indels of length 3\n",
      "44194 indels of length 4\n",
      "18108 indels of length 5\n",
      "9882 indels of length 6\n",
      "4259 indels of length 7\n",
      "2256 indels of length 8\n",
      "1654 indels of length 9\n",
      "866 indels of length 10\n",
      "492 indels of length 11\n",
      "458 indels of length 12\n",
      "206 indels of length 13\n",
      "126 indels of length 14\n",
      "118 indels of length 15\n",
      "75 indels of length 16\n",
      "33 indels of length 17\n",
      "38 indels of length 18\n",
      "17 indels of length 19\n",
      "10 indels of length 20\n",
      "11 indels of length 21\n",
      "7 indels of length 22\n",
      "3 indels of length 23\n",
      "4 indels of length 24\n",
      "7 indels of length 25\n",
      "4 indels of length 26\n",
      "1 indels of length 27\n",
      "1 indels of length 28\n",
      "1 indels of length 30\n",
      "3 indels of length 32\n",
      "2 indels of length 33\n",
      "1 indels of length 34\n",
      "3 indels of length 36\n",
      "1 indels of length 37\n",
      "1 indels of length 43\n",
      "2 indels of length 45\n",
      "1 indels of length 49\n",
      "1 indels of length 52\n",
      "1 indels of length 54\n",
      "1 indels of length 72\n",
      "5 indels of length 113\n",
      "2 indels of length 115\n",
      "1 indels of length 116\n",
      "1 indels of length 148\n",
      "1 indels of length 149\n",
      "1 indels of length 189\n",
      "1 indels of length 236\n",
      "1 indels of length 697\n",
      "1 indels of length 713\n",
      "1 indels of length 763\n",
      "1 indels of length 764\n",
      "1 indels of length 768\n",
      "1 indels of length 1137\n",
      "1 indels of length 1162\n",
      "Indels near homopolymers: 77%\n"
     ]
    }
   ],
   "source": [
    "print(\"========== ONT_2D.fastq.gz  ==========\")\n",
    "\n",
    "genome = next(SeqIO.parse(\"reference.fasta\", \"fasta\")).seq\n",
    "\n",
    "# coverage & errors\n",
    "cov = np.zeros(shape=(len(genome)))\n",
    "SNP = pd.DataFrame(0, index=['A', 'C', 'G', 'T', '-'], columns=['A', 'C', 'G', 'T', '-'])\n",
    "\n",
    "\n",
    "delta = 0\n",
    "indels = Counter() # length -> count\n",
    "homopoly = 0\n",
    "\n",
    "for i, line in enumerate(open(\"ONT.pileup\")):\n",
    "    _, pos, gen, num_aligned, seqs, quals = line.strip().split(\"\\t\")\n",
    "    cov[int(pos)-1] = int(num_aligned)\n",
    "    if int(num_aligned) > 0:\n",
    "        q = 0\n",
    "        while q < len(seqs):\n",
    "            if seqs[q] == '^':\n",
    "                q += 2\n",
    "            elif seqs[q] == '$' or seqs[q] == '*':\n",
    "                q += 1\n",
    "            elif seqs[q] == '.' or seqs[q] == ',':\n",
    "                q += 1\n",
    "            elif seqs[q].upper() in \"AGTC\":\n",
    "                SNP.loc[seqs[q].upper(), gen] += 1\n",
    "                q += 1\n",
    "            elif seqs[q] in \"+-\":\n",
    "                INS = True\n",
    "                if seqs[q] == '-':\n",
    "                    INS = False\n",
    "                num = \"\"\n",
    "                q += 1\n",
    "                while seqs[q] in \"0123456789\":\n",
    "                    num += seqs[q]\n",
    "                    q += 1\n",
    "                num = int(num)\n",
    "                indels[num] += 1\n",
    "                \n",
    "                if gen == genome[int(pos)-2].upper() or gen == genome[int(pos)].upper() or \\\n",
    "                    num > 1 and seqs[q].upper() == seqs[q+1].upper():\n",
    "                    homopoly += num\n",
    "            \n",
    "                \n",
    "                for j in range(num):\n",
    "                    val = seqs[q+j].upper()\n",
    "                    if INS:\n",
    "                        SNP.loc['-', val] += 1\n",
    "                    else:\n",
    "                        SNP.loc[val, '-'] += 1\n",
    "                q += num\n",
    "                \n",
    "                \n",
    "\n",
    "print(\"Number of substitutions:\")\n",
    "print(SNP)\n",
    "total = sum(cov)\n",
    "print(f'mean_read_error={round(SNP.iloc[:4, :].values.sum()/total*100, 5)}%')\n",
    "\n",
    "\n",
    "covs = [np.mean(cov[i:i+1000]) for i in range(0, len(cov), 1000)]\n",
    "plt.figure(figsize=(15, 5))\n",
    "plt.bar(range(len(covs)), height=covs)\n",
    "plt.xticks(range(0, len(covs), 500), range(0, len(cov), 500000))\n",
    "plt.title(\"Mean coverage\")\n",
    "plt.show()\n",
    "print(f'mean_cov={round(np.mean(cov), 2)}, breadth_cov={np.count_nonzero(cov)/len(cov)}')\n",
    "\n",
    "print('Indel stats')\n",
    "for k, v in sorted(indels.items()):\n",
    "    print(f'{v} indels of length {k}')\n",
    "    \n",
    "print(f'Indels near homopolymers: {round(homopoly/sum(indels.values())*100)}%')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "super-immigration",
   "metadata": {},
   "source": [
    "У Oxford Nanopore более длинные инделы, и они гораздо чаще встречаются рядом с гомополимерными участками. У PacBio чаще встречаются mismatches, чем у ONT. Глубина покрытия у ONT больше, чем у PacBio. Процент ошибок на порядок выше, чем в ридах Illumina из дз2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "equal-experiment",
   "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}