checkCompat_110818_1519.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from itertools import groupby\n",
    "from operator import itemgetter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "========== 1\n",
      "should 3\n",
      "writing tmp frameshift\n",
      "========== 68\n",
      "should 6\n",
      "1\n",
      "1 100 68 71 1\n",
      "writing tmp frameshift 6\n",
      "========== 78\n",
      "should 4\n",
      "weird 4\n",
      "========== 200\n",
      "should 2\n",
      "end 2\n",
      "1\n"
     ]
    }
   ],
   "source": [
    "# Now we can proceed to write a comparrison function to parse these exon-intron chains\n",
    "# the function is not expected to be particularly fast\n",
    "# however, what I need is just to get an answer, so...\n",
    "\n",
    "def checkCompat(orf,exons):\n",
    "    # first convert the orf into lists\n",
    "    tmp=orf.split(\"-\")\n",
    "    orfStarts=[int(x) for x in tmp[0].split(\",\")]\n",
    "    orfEnds=[int(x) for x in tmp[-1].split(\",\")]\n",
    "    \n",
    "    # next we need to parse the exons in the same manner\n",
    "    tmp=exons.split(\"-\")\n",
    "    exonStarts=[int(x) for x in tmp[0].split(\",\")]\n",
    "    exonEnds=[int(x) for x in tmp[-1].split(\",\")]\n",
    "    numExons=len(exonStarts)\n",
    "    \n",
    "    # now we need to compare them to one another\n",
    "    # there could be three possibilities which we shall use the following return codes for:\n",
    "    # 0 - ORF fits entirely within the exon\n",
    "    # 1 - section of the ORF is lost without a frameshift\n",
    "    # 2 - an extra section is added to the ORF without causing a framesift\n",
    "    # 3 - there is a frameshift due to an extra section added to the ORF\n",
    "    # 4 - there is a frameshift due to a section removed from the ORF\n",
    "    # 5 - ORF is trimmed at the end - perhaps the exon should be extended???\n",
    "    # 6 - potentially missing start codon - if the beginning of the ORF is trimmed\n",
    "    \n",
    "    assert len(orfStarts)==len(orfEnds),\"wrong CDS chain passed, length don't match: \"+orf\n",
    "    assert len(exonStarts)==len(exonEnds),\"wrong exon chain passed, length don't match: \"+exons\n",
    "#     assert len(orfStarts)<=len(exonStarts),\"There are more segments in ORF than in the exon chain. ORF: \"+orf+\" ; Exons-chain: \"+exons\n",
    "    \n",
    "    frameShift=0\n",
    "    tmpFrameShift=0 # used when novel introns are present within a known ORF - multiple frameShifts need to be evaluated\n",
    "    modification=0\n",
    "    \n",
    "    orfPos=0\n",
    "    exonPos=-1\n",
    "    numOrfSegments=len(orfStarts)\n",
    "    for es,ee in zip(exonStarts,exonEnds):\n",
    "        exonPos+=1\n",
    "        print(\"==========\",es)\n",
    "        curOrfStart=orfStarts[orfPos]\n",
    "        curOrfEnd=orfEnds[orfPos]\n",
    "        \n",
    "        # if this happens before the first ORF - then good\n",
    "        # otherwise we have an insertion and need to check for a frameshift\n",
    "        if ee<curOrfStart:\n",
    "            if orfPos==0:\n",
    "                print(\"continue 0\")\n",
    "                continue\n",
    "            else:\n",
    "                print(\"inserted exon\")\n",
    "                frameShift=((ee+1)-es)%3\n",
    "                if frameShift==0:\n",
    "                    modification=2\n",
    "                else:\n",
    "                    return 3\n",
    "        \n",
    "        # ORF within exon and fits nicely - exon/CDS segment share end\n",
    "        elif es<=curOrfStart and ee==curOrfEnd:\n",
    "            print(\"should 1\")\n",
    "            if orfPos==numOrfSegments-1: # ORF fits and is done\n",
    "                print(\"end 1\")\n",
    "                return modification\n",
    "            else: # ORF continues\n",
    "                if exonPos==numExons-1: # trimmed based on having a shorter exon chain\n",
    "                    return 5\n",
    "                print(\"continue 1\")\n",
    "                orfPos+=1\n",
    "        \n",
    "        # could be the proper end of the ORF or a signal for modification\n",
    "        elif es<=curOrfStart and ee>curOrfEnd:\n",
    "            print(\"should 2\")\n",
    "            if orfPos==numOrfSegments-1: # ORF fits and is done\n",
    "                print(\"end 2\")\n",
    "                return modification\n",
    "            else: # ORF does not match current exon chain\n",
    "                print(\"weird 2\")\n",
    "        \n",
    "        # could be the end of the ORF, or a signal for modification\n",
    "        elif es<=curOrfStart and ee<curOrfEnd:\n",
    "            print(\"should 3\")\n",
    "            # this is another perplexing case\n",
    "            # if an exon exists such that it ends before the end of the current ORF segment\n",
    "            # but there is another exon which fits within the current ORF segment\n",
    "            # then there is a possibility that two frameshifts will cancel out one another\n",
    "            frameShift=(curOrfEnd-ee)%3\n",
    "            if frameShift==0:\n",
    "                print(\"frameshift=0 3\")\n",
    "                modification=1\n",
    "                continue\n",
    "            else:\n",
    "                if exonPos<numExons-1 and exonStarts[exonPos+1]<curOrfEnd: # additional frameshifts need to be evaluated\n",
    "                    print(\"writing tmp frameshift\")\n",
    "                    tmpFrameShift=frameShift\n",
    "                    continue\n",
    "                else:\n",
    "                    return 4\n",
    "        \n",
    "        # frameshift due to exon beginning later\n",
    "        elif es>curOrfStart and ee==curOrfEnd:\n",
    "            print(\"should 4\")\n",
    "            frameShift=((es-curOrfStart)+tmpFrameShift)%3\n",
    "            if orfPos==numOrfSegments-1:\n",
    "                print(\"end 4\")\n",
    "                if frameShift==0:\n",
    "                    return 1\n",
    "                else:\n",
    "                    return 4\n",
    "            else:\n",
    "                print(\"weird 4\")\n",
    "                if frameShift==0:\n",
    "                    orfPos+=1\n",
    "                    modification=1\n",
    "                    continue\n",
    "                else:\n",
    "                    return 4\n",
    "        \n",
    "        # Something\n",
    "        elif es>curOrfStart and ee>curOrfEnd:\n",
    "            print(\"should 5\")\n",
    "            frameShift=(es-curOrfStart)%3\n",
    "            if orfPos==numOrfSegments-1:\n",
    "                print(\"end 5\")\n",
    "                if frameShift==0:\n",
    "                    return 1\n",
    "                else:\n",
    "                    return 4\n",
    "            else:\n",
    "                if frameShift==0:\n",
    "                    print(\"weird 5\")\n",
    "                else:\n",
    "                    return 4 # just abandon - not going to be any worth - it's already a frameshift\n",
    "        \n",
    "        # an intron appeared in the middle of the ORF\n",
    "        elif es>curOrfStart and ee<curOrfEnd:\n",
    "            print(\"should 6\")\n",
    "            # this is another perplexing case\n",
    "            # if an exon exists such that it ends before the end of the current ORF segment\n",
    "            # but there is another exon which fits within the current ORF segment\n",
    "            # then there is a possibility that two frameshifts will cancel out one another\n",
    "\n",
    "            # here frameshift needs to take into account the previously seen tmpFrameShift\n",
    "            print(tmpFrameShift)\n",
    "            frameShift=((curOrfEnd-ee)+((es-curOrfStart)+tmpFrameShift))%3\n",
    "            tmpFrameShift=frameShift\n",
    "            print(curOrfStart,curOrfEnd,es,ee,tmpFrameShift)\n",
    "            if frameShift==0:\n",
    "                print(\"frameshift=0 6\")\n",
    "                modification=1\n",
    "                continue\n",
    "            else:\n",
    "                if exonPos<numExons-1 and exonStarts[exonPos+1]<curOrfEnd: # additional frameshifts need to be evaluated\n",
    "                    print(\"writing tmp frameshift 6\")\n",
    "                    tmpFrameShift=frameShift\n",
    "                    continue\n",
    "                else:\n",
    "                    return 4\n",
    "            \n",
    "        else:\n",
    "            print(\"ELSE:\",es,ee,curOrfStart,curOrfEnd)\n",
    "            break\n",
    "                \n",
    "    return modification\n",
    "    \n",
    "knownExons=\"1,200-100,300\"\n",
    "knownCDS=\"50,200-100,250\"\n",
    "novelExons0=\"1,200,400-100,350,500\" # 0 - the known ORF fits nicely\n",
    "novelExons1=\"1,203-100,300\" # 1 - the known ORF does not fit entirely. However the deletion does not cause a frameshift\n",
    "novelExons2=\"1,150,200-100,161,300\" # 2 - the known ORF does not fit entirely. However, because it is 12 bases that are added - there is no frame shift\n",
    "novelExons3=\"1,150,200-100,162,300\" # 3 - the known ORF does not fit entirely and the insertion causes a frameshift\n",
    "novelExons4=\"1,210-100,300\" # 4 - the known ORF does not fit entirely and the deletion causes a frameshift\n",
    "novelExons5=\"1,200-100,240\"\n",
    "novelExons6=\"60,200-100,300\"\n",
    "novelExons7=\"1-100\" # first segment of the exon chain fits well, but the remaining segments are trimmed\n",
    "\n",
    "\n",
    "# let's now deal with the second batch of test cases\n",
    "novelExons1=\"1,200-97,300\"\n",
    "\n",
    "# test for the case when an extra intron splits an ORF segment in two\n",
    "novelExons1=\"1,78,200-71,100,300\"\n",
    "\n",
    "# to make things even more complicated - what if we introduce two novel introns into this story?\n",
    "# the algorithm should be able to handle them easily\n",
    "knownCDS=\"1,200-100,250\"\n",
    "novelExons1=\"1,68,78,200-63,71,100,300\"\n",
    "\n",
    "# same as before, but a simplified example for ease of manual verification\n",
    "# knownCDS=\"1-15\"\n",
    "# novelExons1=\"1,6,13-3,11,15\"\n",
    "\n",
    "# yet another case when a series of introns interrupt + trimming of the end\n",
    "# novelExons1=\"1,63,78-68,71,99\"\n",
    "\n",
    "# needs to be able to reset the tmpFrameShift\n",
    "\n",
    "# print(checkCompat(knownCDS,knownExons))\n",
    "# print(\"\\n+++++++++++++++\\n\")\n",
    "# print(checkCompat(knownCDS,novelExons0))\n",
    "# print(\"\\n+++++++++++++++\\n\")\n",
    "print(checkCompat(knownCDS,novelExons1))\n",
    "# print(\"\\n+++++++++++++++\\n\")\n",
    "# print(checkCompat(knownCDS,novelExons2))\n",
    "# print(\"\\n+++++++++++++++\\n\")\n",
    "# print(checkCompat(knownCDS,novelExons3))\n",
    "# print(\"\\n+++++++++++++++\\n\")\n",
    "# print(checkCompat(knownCDS,novelExons4))\n",
    "# print(\"\\n+++++++++++++++\\n\")\n",
    "# print(checkCompat(knownCDS,novelExons5))\n",
    "# print(\"\\n+++++++++++++++\\n\")\n",
    "# print(checkCompat(knownCDS,novelExons6))\n",
    "# print(\"\\n+++++++++++++++\\n\")\n",
    "# print(checkCompat(knownCDS,novelExons7))\n",
    "\n",
    "# needs to also account for trimming, when the ORF has more segments than the exon-chain to which it is compared\n",
    "\n",
    "# So far we've tested cases when it's the exonStart that affects the match\n",
    "# now need to test the reverse - exonEnd doesn't fit\n",
    "# tests should be performed for each of the 6 cases as well\n",
    "\n",
    "# next we also need to test for the presence of multiple modifications and whether they are reported as needed\n",
    "\n",
    "# also need to test longer exon chains\n",
    "\n",
    "# also need to test longer CDS chains"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def checkCompat(orf,exons):\n",
    "    # first convert the orf into lists\n",
    "    tmp=orf.split(\"-\")\n",
    "    orfStarts=[int(x) for x in tmp[0].split(\",\")]\n",
    "    orfEnds=[int(x) for x in tmp[-1].split(\",\")]\n",
    "    \n",
    "    # next we need to parse the exons in the same manner\n",
    "    tmp=exons.split(\"-\")\n",
    "    exonStarts=[int(x) for x in tmp[0].split(\",\")]\n",
    "    exonEnds=[int(x) for x in tmp[-1].split(\",\")]\n",
    "    \n",
    "    numExons=len(exonStarts)-1\n",
    "    numOrfSegments=len(orfStarts)-1\n",
    "    \n",
    "    assert len(orfStarts)==len(orfEnds),\"wrong CDS chain passed, length don't match: \"+orf\n",
    "    assert len(exonStarts)==len(exonEnds),\"wrong exon chain passed, length don't match: \"+exons\n",
    "    \n",
    "    orfPos=0\n",
    "    exonPos=0\n",
    "    \n",
    "    code=\"\"\n",
    "    orfNotStarted=True # orf has not been reached yet\n",
    "    \n",
    "    while True:\n",
    "        es=exonStarts[exonPos]\n",
    "        ee=exonEnds[exonPos]\n",
    "        os=orfStarts[orfPos]\n",
    "        oe=orfEnds[orfPos]\n",
    "        # skip to the exon that overlaps the start of the ORF if such exists\n",
    "        if ee<os:\n",
    "            if exonPos==numExons and orfNotStarted:\n",
    "                return \"9999999\"\n",
    "            exonPos+=1\n",
    "        else:\n",
    "            print(\"es\",es,ee,os,oe)\n",
    "            if es>=os:\n",
    "                print(\"if 1\")\n",
    "                if ee==oe:\n",
    "                    print(\"if 2\")\n",
    "                    if orfPos==numOrfSegments:\n",
    "                        return code+\";0\"\n",
    "                    exonPos+=1\n",
    "                    orfPos+=1\n",
    "                    continue\n",
    "                elif ee>oe:\n",
    "                    print(\"elif 2\")\n",
    "                    if orfPos==numOrfSegments:\n",
    "                        return code\n",
    "                    else: # exon extends past the end of the orf segment - orf3-chain1\n",
    "                        # deal with the frameshift here\n",
    "                        \n",
    "                        # if exon fills the entire ORF intron -> os(next orf)-oe(current orf) -> change exon start to the end of orf intron (orf3-chain1)\n",
    "                        if ee>orfStarts[orfPos+1]:\n",
    "                            print(\"if 3\")\n",
    "                            numAdded=(orfStarts[orfPos+1]-1)-oe\n",
    "                            frameshift=numAdded%3\n",
    "                            if frameshift==0:\n",
    "                                code+=\";2\"\n",
    "                            else:\n",
    "                                code+=\";3\"\n",
    "                            \n",
    "                            exonStarts[exonPos]=orfStarts[orfPos+1]\n",
    "                            orfPos+=1\n",
    "                        # else if exon fills part of the ORF intron -> ee-oe -> (orf3-chain2)\n",
    "                        else:\n",
    "                            print(\"else 3\")\n",
    "                            numAdded=ee-oe\n",
    "                            frameshift=numAdded%3\n",
    "                            if frameshift==0:\n",
    "                                code+=\";2\"\n",
    "                            else:\n",
    "                                code+=\";3\"\n",
    "                            \n",
    "                            exonPos+=1\n",
    "                            orfPos+=1\n",
    "                else:\n",
    "                    print(\"else 4\")\n",
    "                    if exonPos==numExons:\n",
    "                        return code+\";5\"\n",
    "                    else:\n",
    "                        print(\"else 5\")\n",
    "                        # if orf fills the entire space between two exons\n",
    "                        if oe>exonStarts[exonPos+1]:\n",
    "                            print(\"if 6\")\n",
    "                            numRemoved=(exonStarts[exonPos+1]-1)-ee\n",
    "                            frameshift=numRemoved%3\n",
    "                            if frameshift==0:\n",
    "                                code+=\";1\"\n",
    "                            else:\n",
    "                                code+=\";4\"\n",
    "                            \n",
    "                            orfStarts[orfPos]=exonStarts[exonPos]\n",
    "                            exonPos+=1\n",
    "                        else:\n",
    "                            print(\"else 6\")\n",
    "                            numRemoved=oe-ee\n",
    "                            frameshift=numremoved%3\n",
    "                            if frameshift==0:\n",
    "                                code+=\";1\"\n",
    "                            else:\n",
    "                                code+=\";4\"\n",
    "                            \n",
    "                            exonPos+=1\n",
    "                            orfPos+=1\n",
    "\n",
    "            else:\n",
    "                return \"adsjflkads\"\n",
    "                "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 160,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(1, 2), (8, 9), (15, 15)]"
      ]
     },
     "execution_count": 160,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# from https://codereview.stackexchange.com/questions/178427/given-2-disjoint-sets-of-intervals-find-the-intersections\n",
    "\n",
    "def getIntersection(interval_1, interval_2):\n",
    "    start = max(interval_1[0], interval_2[0])\n",
    "    end = min(interval_1[1], interval_2[1])\n",
    "    if start <= end:\n",
    "        return (start, end)\n",
    "    return None\n",
    "\n",
    "def getIntersection2(interval_1,interval_2):\n",
    "    left_start=\n",
    "\n",
    "def return_intersections(intervals1, intervals2):\n",
    "    iter1 = iter(intervals1)\n",
    "    iter2 = iter(intervals2)\n",
    "\n",
    "    interval1 = next(iter1)\n",
    "    interval2 = next(iter2)\n",
    "\n",
    "    while True:\n",
    "        intersection = getIntersection(interval1, interval2)\n",
    "        if intersection:\n",
    "            yield intersection\n",
    "            try:\n",
    "                if intersection[1] == interval1[1]:\n",
    "                    interval1 = next(iter1)\n",
    "                else:\n",
    "                    interval2 = next(iter2)\n",
    "            except StopIteration:\n",
    "                return\n",
    "        \n",
    "        else: # need to handle this last case so it does not go into an infinite loop\n",
    "            return\n",
    "        \n",
    "        try:\n",
    "            while interval1[0] > interval2[1]:\n",
    "                interval2 = next(iter2)\n",
    "            while interval2[0] > interval1[1]:\n",
    "                interval1 = next(iter1)\n",
    "        except StopIteration:\n",
    "            return\n",
    "        \n",
    "orf1=[(1,18)]\n",
    "orf3=[(1,3),(7,15)]\n",
    "chain1=[(1,18)]\n",
    "chain2=[(1,6),(10,18)]\n",
    "chain17=[(1,2),(8,9),(15,18)]\n",
    "list(return_intersections(orf3,chain17))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{1, 2, 3, 7, 8, 9, 10, 11, 12, 13, 14, 15}\n",
      "{1, 2, 9, 15}\n",
      "{3, 7, 8, 10, 11, 12, 13, 14}\n",
      "set()\n",
      "[[3], [7, 8], [10, 11, 12, 13, 14]]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "2"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# what if we try sets\n",
    "# resChain=[]\n",
    "# for t in chain2:\n",
    "#     resChain=resChain+list(range(t[0],t[1]))\n",
    "# resOrf=[]\n",
    "# for t in orf3:\n",
    "#     resOrf=resOrf+list(range(t[0],t[1]))\n",
    "    \n",
    "# print(set(resOrf)-set(resChain))\n",
    "# print(set(resChain)-set(resOrf))\n",
    "\n",
    "def checkCompat(orf,chain):\n",
    "    tmp=orf.split(\"-\")\n",
    "    orfStarts=[int(x) for x in tmp[0].split(\",\")]\n",
    "    orfEnds=[int(x) for x in tmp[-1].split(\",\")]\n",
    "    \n",
    "    # next we need to parse the exons in the same manner\n",
    "    tmp=chain.split(\"-\")\n",
    "    chainStarts=[int(x) for x in tmp[0].split(\",\")]\n",
    "    chainEnds=[int(x) for x in tmp[-1].split(\",\")]\n",
    "    \n",
    "    orf=[]\n",
    "    for t in zip(orfStarts,orfEnds):\n",
    "        orf+=list(range(t[0],t[1]+1))\n",
    "    \n",
    "    chain=[]\n",
    "    for t in zip(chainStarts,chainEnds):\n",
    "        chain+=list(range(t[0],t[1]+1))\n",
    "    \n",
    "    # now remove anything in the chain that is before or after the start and end of the ORF\n",
    "    chain=[x for x in chain if x>=orf[0] and x<=orf[-1]]\n",
    "    \n",
    "    orf=set(orf)\n",
    "    print(orf)\n",
    "    chain=set(chain)\n",
    "    print(chain)\n",
    "    \n",
    "    # now get set differences\n",
    "    orf_chain=orf-chain\n",
    "    chain_orf=chain-orf\n",
    "    print(orf_chain)\n",
    "    print(chain_orf)\n",
    "    \n",
    "    if len(orf_chain)==0 and len(chain_orf)==0:\n",
    "        print(\"hooray\")\n",
    "        return \"0\"\n",
    "    \n",
    "    allMis=[]\n",
    "    for k, g in groupby(enumerate(orf_chain), lambda ix : ix[0] - ix[1]):\n",
    "        allMis.append(list(map(itemgetter(1), g)))\n",
    "    for k, g in groupby(enumerate(chain_orf), lambda ix : ix[0] - ix[1]):\n",
    "        allMis.append(list(map(itemgetter(1), g)))\n",
    "        \n",
    "    print(allMis)\n",
    "    # now just need to figure out what kind of frameshifts we've got here\n",
    "    offFrame=0 # number of extra or missing bases that do result in a frameshift\n",
    "    onFrame=0 # number of extra or missing bases that do not result in a frameshift\n",
    "    \n",
    "    # for now this will be a very easy thing to do\n",
    "    # would be better to have a method for computing how many bases are in off-frame vs on-frame regions\n",
    "    for mis in allMis:\n",
    "        if not len(mis)%3==0:\n",
    "            return 2\n",
    "    return 1\n",
    "    \n",
    "orf3=\"1,7-3,15\"\n",
    "chain17=\"1,9,15-2,9,18\"\n",
    "checkCompat(orf3,chain17)\n",
    "\n",
    "# need to figure out a way to replace the codes with percentage as follows:\n",
    "# (number of bases in off-frame)/(total number of bases in the ORF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{1, 2, 3, 7, 8, 9, 10, 11, 12, 13, 14, 15}\n",
      "{1, 2, 9, 15}\n",
      "{3, 7, 8, 10, 11, 12, 13, 14}\n",
      "set()\n",
      "[[3], [7, 8], [10, 11, 12, 13, 14]]\n",
      "[]\n",
      "[[], [], []]\n",
      "{1, 2, 3, 7, 8, 9, 10, 11, 12, 13, 14, 15}\n",
      "{1, 2, 9, 15}\n"
     ]
    }
   ],
   "source": [
    "# what if we try sets\n",
    "# resChain=[]\n",
    "# for t in chain2:\n",
    "#     resChain=resChain+list(range(t[0],t[1]))\n",
    "# resOrf=[]\n",
    "# for t in orf3:\n",
    "#     resOrf=resOrf+list(range(t[0],t[1]))\n",
    "    \n",
    "# print(set(resOrf)-set(resChain))\n",
    "# print(set(resChain)-set(resOrf))\n",
    "\n",
    "def checkCompat(orf,chain):\n",
    "    tmp=orf.split(\"-\")\n",
    "    orfStarts=[int(x) for x in tmp[0].split(\",\")]\n",
    "    orfEnds=[int(x) for x in tmp[-1].split(\",\")]\n",
    "    \n",
    "    # next we need to parse the exons in the same manner\n",
    "    tmp=chain.split(\"-\")\n",
    "    chainStarts=[int(x) for x in tmp[0].split(\",\")]\n",
    "    chainEnds=[int(x) for x in tmp[-1].split(\",\")]\n",
    "    \n",
    "    orf=[]\n",
    "    for t in zip(orfStarts,orfEnds):\n",
    "        orf+=list(range(t[0],t[1]+1))\n",
    "    \n",
    "    chain=[]\n",
    "    for t in zip(chainStarts,chainEnds):\n",
    "        chain+=list(range(t[0],t[1]+1))\n",
    "    \n",
    "    # now remove anything in the chain that is before or after the start and end of the ORF\n",
    "    chain=[x for x in chain if x>=orf[0] and x<=orf[-1]]\n",
    "    \n",
    "    orf=set(orf)\n",
    "    print(orf)\n",
    "    chain=set(chain)\n",
    "    print(chain)\n",
    "    \n",
    "    # now get set differences\n",
    "    orf_chain=orf-chain\n",
    "    chain_orf=chain-orf\n",
    "    print(orf_chain)\n",
    "    print(chain_orf)\n",
    "    \n",
    "    if len(orf_chain)==0 and len(chain_orf)==0:\n",
    "        print(\"hooray\")\n",
    "        return \"0\"\n",
    "    \n",
    "    allMis=[]\n",
    "    allMis_orf_chain=[]\n",
    "    allMis_chain_orf=[]\n",
    "    for k, g in groupby(enumerate(orf_chain), lambda ix : ix[0] - ix[1]):\n",
    "        allMis_orf_chain.append(list(map(itemgetter(1), g)))\n",
    "        allMis.append(list(map(itemgetter(1),g)))\n",
    "    for k, g in groupby(enumerate(chain_orf), lambda ix : ix[0] - ix[1]):\n",
    "        allMis_chain_orf.append(list(map(itemgetter(1), g)))\n",
    "        allMis.append(list(map(itemgetter(1),g)))\n",
    "        \n",
    "    print(allMis_orf_chain)\n",
    "    print(allMis_chain_orf)\n",
    "    print(allMis)\n",
    "    # now just need to figure out what kind of frameshifts we've got here\n",
    "    \n",
    "    # let's try to make it better and compute percentage of off-frame bases instead\n",
    "    \n",
    "    # one of the first differences should be:\n",
    "    # we need to compute what's missing from the ORF\n",
    "    # and what is added to the ORF separately\n",
    "    # meaning two different lists - one as orf.difference(chain)\n",
    "    # and the other as chain.difference(orf)\n",
    "    \n",
    "    # If a frameshift does happen - we want to know how many bases of the ORF are affected\n",
    "    # basically, what we need to do instead of computing frameshifts\n",
    "    # we want to know how many bases are different due to the given modifications due to the peptide sequences\n",
    "    \n",
    "    # one thing to take into account is that the two sets are guaranteed to be non-overlapping\n",
    "    # we shall take advantage of that, since we know we can order them\n",
    "    \n",
    "    # first thing first - figure out whether the ORF start is ocvered by the exon:    \n",
    "    if not next(iter(orf))==next(iter(chain)): # if exon does not cover the start codon\n",
    "        return 555555\n",
    "    \n",
    "    offFrame=0 # number of extra or missing bases that do result in a frameshift\n",
    "    onFrame=0 # number of extra or missing bases that do not result in a frameshift\n",
    "    \n",
    "    removed=False # whether the elements in the current array have been removed or added to the original ORF\n",
    "    curChain=[]\n",
    "    while True: # parsing through all modifications\n",
    "        # first figure out which one should be evaluated now: allMis_orf_chain or the allMis_chain_orf\n",
    "        if not allMis_orf_chain and not allMis_chain_orf:\n",
    "            return\n",
    "        elif not allMis_orf_chain: # allMis_orf_chain is empty\n",
    "            curChain=allMis_chain_orf.pop(0)\n",
    "            removed=False\n",
    "        elif not allMis_chain_orf:# allMis_chain_orf is empty\n",
    "            curChain=allMis_orf_chain.pop(0)\n",
    "            removed=True\n",
    "        elif allMis_chain_orf[0]<allMis_orf_chain[0]: # present in the new transcript but not in the orf\n",
    "            curChain=allMis_chain_orf.pop(0)\n",
    "            removed=False\n",
    "        else: # present in the ORF but not the exon chain\n",
    "            curChain=allMis_orf_chain.pop(0)\n",
    "            removed=True\n",
    "        \n",
    "        # now to proceed with the current chain\n",
    "        # in order to evaluate what overlapping sections of the ORF are being \n",
    "        \n",
    "orf3=\"1,7-3,15\"\n",
    "chain17=\"1,9,15-2,9,18\"\n",
    "checkCompat(orf3,chain17)\n",
    "\n",
    "# need to figure out a way to replace the codes with percentage as follows:\n",
    "# (number of bases in off-frame)/(total number of bases in the ORF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "es 1 6 1 3\n",
      "if 1\n",
      "elif 2\n",
      "else 3\n",
      "es 10 18 7 15\n",
      "if 1\n",
      "elif 2\n",
      "2;1 : ;2\n"
     ]
    }
   ],
   "source": [
    "# orf3:    1+  +  +3           7+  +  +   +  +  +   +  +  +15\n",
    "# orf2:              4+  +  +   +  +  +   +  +  +   +  +  +15\n",
    "# orf1:    1+  +  +   +  +  +   +  +  +   +  +  +   +  +  +   +  +  +18\n",
    "# exons:   |-  -  -  |-  -  -  |-  -  -  |-  -  -  |-  -  -  |-  -  -  |\n",
    "# pos:      1  2  3   4  5  6   7  8  9   10 11 12  13 14 15  16 17 18\n",
    "# chain1:  1=  =  =   =  =  =   =  =  =   =  =  =   =  =  =   =  =  =18\n",
    "# chain2:  1=  =  =   =  =  =6          10=  =  =   =  =  =   =  =  =18 # compatible - deleted ORF segment\n",
    "\n",
    "orf3=\"1,7-3,15\"\n",
    "orf2=\"4-15\"\n",
    "orf1=\"1-18\"\n",
    "\n",
    "chain1=\"1-18\"\n",
    "chain2=\"1,10-6,18\"\n",
    "\n",
    "# print(\"0 :\",checkCompat(orf1,chain1))\n",
    "# print(\"1 :\",checkCompat(orf1,chain2))\n",
    "# print(\"0 :\",checkCompat(orf2,chain1))\n",
    "# print(\"1 :\",checkCompat(orf2,chain2))\n",
    "# print(\"2 :\",checkCompat(orf3,chain1))\n",
    "print(\"2;1 :\",checkCompat(orf3,chain2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "metadata": {},
   "outputs": [],
   "source": [
    "# orf3:    1+  +  +3           7+  +  +   +  +  +   +  +  +15\n",
    "# orf2:              4+  +  +   +  +  +   +  +  +   +  +  +15\n",
    "# orf1:    1+  +  +   +  +  +   +  +  +   +  +  +   +  +  +   +  +  +18\n",
    "# exons:   |-  -  -  |-  -  -  |-  -  -  |-  -  -  |-  -  -  |-  -  -  |\n",
    "# pos:      1  2  3   4  5  6   7  8  9   10 11 12  13 14 15  16 17 18\n",
    "# chain3:  1=  =  =   =  =5             10=  =  =   =  =  =   =  =  =18 # frameshift(1) - deleted ORF segment\n",
    "# chain4:  1=  =  =   =  =5                11=  =   =  =  =   =  =  =18 # frameshift(2) - deleted ORF segment\n",
    "# chain5:  1=  =  =   =  =5                   12=   =  =  =   =  =  =18 # compatible - deleted ORF segment\n",
    "\n",
    "orf3=\"1,7-3,15\"\n",
    "orf2=\"4-15\"\n",
    "orf1=\"1-18\"\n",
    "\n",
    "chain3=\"1,10-5,18\"\n",
    "chain4=\"1,11-5,18\"\n",
    "chain5=\"1,12-5,18\"\n",
    "\n",
    "# orf3:    1+  +  +3           7+  +  +   +  +  +   +  +  +15\n",
    "# orf2:              4+  +  +   +  +  +   +  +  +   +  +  +15\n",
    "# orf1:    1+  +  +   +  +  +   +  +  +   +  +  +   +  +  +   +  +  +18\n",
    "# exons:   |-  -  -  |-  -  -  |-  -  -  |-  -  -  |-  -  -  |-  -  -  |\n",
    "# pos:      1  2  3   4  5  6   7  8  9   10 11 12  13 14 15  16 17 18\n",
    "# chain6:  1=  =  =   =  =  =6             11=  =   =  =  =   =  =  =18 # frameshift(1) - deleted ORF segment\n",
    "# chain7:  1=  =  =   =  =5                11=  =   =  =  =   =  =  =18 # frameshift(2) - deleted ORF segment\n",
    "# chain8:  1=  =  =   =4                   11=  =   =  =  =   =  =  =18 # compatible - deleted ORF segment\n",
    "\n",
    "orf3=\"1,7-3,15\"\n",
    "orf2=\"4-15\"\n",
    "orf1=\"1-18\"\n",
    "\n",
    "chain6=\"1,11-6,18\"\n",
    "chain7=\"1,11-5,18\"\n",
    "chain8=\"1,11-4,18\"\n",
    "\n",
    "# orf3:    1+  +  +3           7+  +  +   +  +  +   +  +  +15\n",
    "# orf2:              4+  +  +   +  +  +   +  +  +   +  +  +15\n",
    "# orf1:    1+  +  +   +  +  +   +  +  +   +  +  +   +  +  +   +  +  +18\n",
    "# exons:   |-  -  -  |-  -  -  |-  -  -  |-  -  -  |-  -  -  |-  -  -  |\n",
    "# pos:      1  2  3   4  5  6   7  8  9   10 11 12  13 14 15  16 17 18\n",
    "# chain9:  1=  =  =3           7=  =  =9          13=  =  =   =  =  =18 # compatible - deleted ORF segment\n",
    "# chain10: 1=  =  =3              8=  =9          13=  =  =   =  =  =18 # frameshift(1) - deleted ORF segment\n",
    "# chain11: 1=  =  =3                 9=9          13=  =  =   =  =  =18 # frameshift(2) - deleted ORF segment\n",
    "# chain12: 1=  =2                    9=9          13=  =  =   =  =  =18 # compatible - deleted ORF segment\n",
    "# chain13: 1=  =  =3              8=  =9             14=  =   =  =  =18 # frameshift(2) - deleted ORF segment\n",
    "# chain14: 1=  =  =3                 9=9                15=   =  =  =18 # frameshift(1) - deleted ORF segment\n",
    "# chain15: 1=  =2 =3              8=  =9                15=   =  =  =18 # compatible - deleted ORF segment\n",
    "# chain16: 1=  =2                 8=  =9                15=   =  =  =18 # incompatible(1) - deleted ORF segment\n",
    "# chain17: 1=1                       9=9                15=   =  =  =18 # compatible - deleted ORF segment\n",
    "# chain18: 1=  =2              7=  =  =9                15=   =  =  =18 # compatible - deleted ORF segment\n",
    "# chain19: 1=  =2              7=  =  =9             14=  =   =  =  =18 # compatible - deleted ORF segment\n",
    "\n",
    "orf3=\"1,7-3,15\"\n",
    "orf2=\"4-15\"\n",
    "orf1=\"1-18\"\n",
    "\n",
    "chain9=\"1,7,13-3,9,18\"\n",
    "chain10=\"1,8,13-3,9,18\"\n",
    "chain11=\"1,9,13-3,9,18\"\n",
    "chain12=\"1,9,13-2,9,18\"\n",
    "chain13=\"1,8,14-3,9,18\"\n",
    "chain14=\"1,9,15-3,9,18\"\n",
    "chain15=\"1,8,15-3,9,18\"\n",
    "chain16=\"1,8,15-2,9,18\"\n",
    "chain17=\"1,9,15-1,9,18\"\n",
    "chain18=\"1,7,15-2,9,18\"\n",
    "chain19=\"1,7,14-2,9,18\""
   ]
  }
 ],
 "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
}