demo.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "def merge_overlap(intervals):\n",
    "    sorted_by_lower_bound = sorted(intervals, key=lambda tup: tup[0])\n",
    "    merged = []\n",
    "    for higher in sorted_by_lower_bound:\n",
    "        if not merged:\n",
    "            merged.append(higher)\n",
    "        else:\n",
    "            lower = merged[-1]\n",
    "            # test for intersection between lower and higher:\n",
    "            # we know via sorting that lower[0] <= higher[0]\n",
    "            if higher[0] <= lower[1]:\n",
    "                upper_bound = max(lower[1], higher[1])\n",
    "                merged[-1] = (lower[0], upper_bound)  # replace by merged interval\n",
    "            else:\n",
    "                merged.append(higher)\n",
    "    return merged"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "df_ann = pd.read_csv('/Users/juan/Documents/manu/dev/vms/bio/data/wheatmitedemo.gff', index_col=False, sep='\\t', header=None, comment='#')\n",
    "df_ann.columns = ['seqname', 'source', 'feature', 'start', 'end', 'score', 'strand', 'frame', 'attribute']\n",
    "\n",
    "df_ref = pd.read_csv('/Users/juan/Documents/manu/dev/vms/bio/MITETracker/extras/data/iwgsc.gff3', index_col=False, sep='\\t', header=None, comment='#')\n",
    "df_ref.columns = ['seqname', 'source', 'feature', 'start', 'end', 'score', 'strand', 'frame', 'attribute']\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "total = {}\n",
    "for k, v in df_ref.iterrows():\n",
    "    total[v.seqname] = v.end\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "overlapped = {}\n",
    "for k, v in df_ann.iterrows():\n",
    "    overlapped.setdefault(v.seqname, []).append((v.start, v.end))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "ename": "KeyError",
     "evalue": "'chrUn'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyError\u001b[0m                                  Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-31-b65886887649>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      2\u001b[0m \u001b[0msequences\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdf_ref\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mseqname\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0munique\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtolist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mseqname\u001b[0m \u001b[0;32min\u001b[0m \u001b[0msequences\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m     \u001b[0moverlapped\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mseqname\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmerge_overlap\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0moverlapped\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mseqname\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mKeyError\u001b[0m: 'chrUn'"
     ]
    }
   ],
   "source": [
    "\n",
    "sequences = df_ref.seqname.unique().tolist()\n",
    "for seqname in sequences:\n",
    "    overlapped[seqname] = merge_overlap(overlapped[seqname])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "seqname chr1A len: 594102056 covered: 12699 percentage: 0.0021375115389265713 \n",
      "seqname chr1B len: 689851870 covered: 10545 percentage: 0.0015285890288302038 \n",
      "seqname chr1D len: 495453186 covered: 7817 percentage: 0.0015777474483734575 \n",
      "seqname chr2A len: 780798557 covered: 9181 percentage: 0.0011758474599742376 \n",
      "seqname chr2B len: 801256715 covered: 14201 percentage: 0.0017723408408502387 \n",
      "seqname chr2D len: 651852609 covered: 13258 percentage: 0.0020338953648339235 \n",
      "seqname chr3A len: 750843639 covered: 11022 percentage: 0.0014679487748846735 \n",
      "seqname chr3B len: 830829764 covered: 9428 percentage: 0.0011347691679471415 \n",
      "seqname chr3D len: 615552423 covered: 10318 percentage: 0.0016762179165364117 \n",
      "seqname chr4A len: 744588157 covered: 10296 percentage: 0.0013827778353987439 \n",
      "seqname chr4B len: 673617499 covered: 13441 percentage: 0.001995346026484386 \n",
      "seqname chr4D len: 509857067 covered: 5480 percentage: 0.001074811031304191 \n",
      "seqname chr5A len: 709773743 covered: 14118 percentage: 0.0019890845694470836 \n",
      "seqname chr5B len: 713149757 covered: 10177 percentage: 0.001427049494177981 \n",
      "seqname chr5D len: 566080677 covered: 9712 percentage: 0.0017156565123313686 \n",
      "seqname chr6A len: 618079260 covered: 7322 percentage: 0.00118463771135113 \n",
      "seqname chr6B len: 720988478 covered: 8918 percentage: 0.0012369129704732952 \n",
      "seqname chr6D len: 473592718 covered: 6949 percentage: 0.0014672945203519788 \n",
      "seqname chr7A len: 736706236 covered: 12113 percentage: 0.0016442103253758802 \n",
      "seqname chr7B len: 750620385 covered: 21142 percentage: 0.002816603495254129 \n",
      "seqname chr7D len: 638686055 covered: 12063 percentage: 0.0018887213687482185 \n"
     ]
    },
    {
     "ename": "KeyError",
     "evalue": "'chrUn'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyError\u001b[0m                                  Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-32-f8e26aabad77>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      2\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mseqname\u001b[0m \u001b[0;32min\u001b[0m \u001b[0msequences\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 3\u001b[0;31m     \u001b[0mcurrent\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0moverlapped\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mseqname\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      4\u001b[0m     \u001b[0mtotal_length\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m     \u001b[0;32mfor\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mstart\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mend\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mcurrent\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyError\u001b[0m: 'chrUn'"
     ]
    }
   ],
   "source": [
    "\n",
    "for seqname in sequences:\n",
    "    current = overlapped[seqname]\n",
    "    total_length = 0\n",
    "    for (start, end) in current:\n",
    "        total_length += abs(end - start)\n",
    "    current_ref = df_ref[(df_ref.seqname == seqname)].iloc[0]\n",
    "    ref_len = abs(current_ref.end - current_ref.start)\n",
    "    perc = total_length * 100 / ref_len\n",
    "    print('seqname %s len: %s covered: %s percentage: %s ' % (seqname, ref_len, total_length, perc,))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": 2
}