syntenicCoordinateShift.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Syntenic Coordinate Shift"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Project information"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "This notebook is essentially a set of functions that will allow for the conversion of one genomes coordinates to another via their syntenic relationship (if it exists)\n",
    "\n",
    "The program utilizes the syntenic alignment outputs from the program Satsuma to shift coordinates from one genome to another.  \n",
    "Syntenic blocks between the two input genomes are found and defined by Satsuma and written out to two main files...  \n",
    "**satsuma_summary.chained.out** - The most succinct version of the sytenic blocks with only the coordinates and orientation of the two regions with respect to eachother.  \n",
    "**MergeXCorrMatches.chained.out** - The same information as the summary.chained.out file, but contains the actual alignment that was found between the two regions. **This is the file that is used as input into the coordinate shifting algorithm.**\n",
    "\n",
    "The program works in two steps (with a bit more detailed descriptions later in the notebook) -  \n",
    "1. Read the reference genome scaffold/chromosome sequences into memory (fastaFilePaths)  \n",
    "   Read the MergeXCorrMatches.chained.out file into memory in the form of two interval trees (one for each target and queary genomes) with adjusted genomic coordinates based on the reference genomes  \n",
    "2. Queary any given Chromosome,position in one genome to find the sytenic region and corresponding base pair in the other genome. There are two functions that accomplish this...  \n",
    "-The first is **quearyTree** wich will identify the sytenic region/block(s) that the position falls within.  \n",
    "-The second is **quearySynBlock** which will take the SynBlock object found in the first function, and find the coordinate/base pair information found between genome1 and genome2.\n",
    "\n",
    "Essentially, the first step of the algorithm takes about 5-10 minutes to run depending on the machine but only needs to be run once eachtime the notebook/script is opened. Then all point quearies are super fast.  \n",
    "The alternative to this would be a database, but there are numerous reasons why we chose not to go that route"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-10-21T13:54:29.201085Z",
     "start_time": "2020-10-21T13:54:28.895967Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": [
    "import time\n",
    "import random\n",
    "import os\n",
    "import intervaltree\n",
    "import pickle"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Satsuma homology to memory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-10-21T13:54:31.296325Z",
     "start_time": "2020-10-21T13:54:31.248309Z"
    },
    "code_folding": [
     0,
     17,
     42,
     59,
     70,
     77,
     112,
     189
    ],
    "hidden": true
   },
   "outputs": [],
   "source": [
    "def readFastaToMem(fastaFile):\n",
    "    inFile = open(fastaFile,'r')\n",
    "    chromDict = {}\n",
    "    for line in inFile:\n",
    "        line = line.strip('\\r').strip('\\n')\n",
    "        if line[0] == \">\":\n",
    "            chrom = line[1:]\n",
    "            chromDict.update({ chrom:[] })\n",
    "        else:\n",
    "            chromDict[chrom].append(line)\n",
    "    ##############\n",
    "    inFile.close()\n",
    "    for c in chromDict:\n",
    "        chromDict[c] = ''.join(chromDict[c])\n",
    "    ################\n",
    "    return chromDict\n",
    "\n",
    "def parseIntervalLines(inFile,line1,line2):\n",
    "    cols1 = line1.strip('\\r').strip('\\n').split(' ')\n",
    "    cols2 = line2.strip('\\r').strip('\\n').split(' ')\n",
    "    \n",
    "    orientation = cols1[7]\n",
    "    length = int(cols1[9])\n",
    "    quearyChrom,targetChrom = cols1[1],cols1[5]\n",
    "    quearyInterval = cols1[2].strip(']').strip('[').split('-')\n",
    "    quearyStart,quearyStop = int(quearyInterval[0]),int(quearyInterval[1])\n",
    "    targetInterval = cols1[6].strip(']').strip('[').split('-')\n",
    "    targetStart,targetStop = int(targetInterval[0]),int(targetInterval[1])\n",
    "    pID = float(cols2[4])\n",
    "    #if abs(targetStart-targetStop) != abs(quearyStart-quearyStop):\n",
    "    #    print(line)\n",
    "\n",
    "    \n",
    "    synBlockObj = SynBlock(quearyChrom,quearyStart,quearyStop,\n",
    "                           targetChrom,targetStart,targetStop,\n",
    "                           orientation,length,pID)\n",
    "    \n",
    "    inFile.readline()\n",
    "    inFile.readline()\n",
    "    \n",
    "    return inFile,synBlockObj\n",
    "\n",
    "class SynBlock:\n",
    "    def __init__(self,quearyName,quearyBegin,quearyEnd,targetName,targetBegin,targetEnd,orientation,length,pID):\n",
    "        self.quearyName = quearyName\n",
    "        self.quearyBegin = quearyBegin\n",
    "        self.quearyEnd = quearyEnd\n",
    "        self.targetName = targetName\n",
    "        self.targetBegin = targetBegin\n",
    "        self.targetEnd = targetEnd\n",
    "        self.orientation = orientation\n",
    "        self.length = length\n",
    "        self.pID = pID\n",
    "    ############################\n",
    "    def getSequences(self,seqs):\n",
    "        self.quearySeq = ''.join(seqs[0])\n",
    "        self.linkSeq = ''.join(seqs[1])\n",
    "        self.targetSeq = ''.join(seqs[2])\n",
    "\n",
    "def processAlignment(inFile,alignment):\n",
    "    alignment[1].append(inFile.readline().strip('\\r').strip('\\n'))\n",
    "    alignment[2].append(inFile.readline().strip('\\r').strip('\\n'))\n",
    "    \n",
    "    inFile.readline()\n",
    "    inFile.readline()\n",
    "    inFile.readline()\n",
    "    inFile.readline()\n",
    "    \n",
    "    return inFile,alignment\n",
    "\n",
    "def revCompSatsuma(seq):\n",
    "    compDict = { \"A\":\"T\",\"T\":\"A\",\"G\":\"C\",\"C\":\"G\",\n",
    "                 \"a\":\"t\",\"t\":\"a\",\"g\":\"c\",\"c\":\"g\",\n",
    "                 \"N\":\"N\",\"n\":\"n\",\"-\":\"-\" }\n",
    "    ##################################################\n",
    "    return ''.join([ compDict[s] for s in seq ][::-1])\n",
    "\n",
    "def preProcessIntervalCoords(synBlockObj,quearyRefDict,targetRefDict,howFar=2000):\n",
    "    qSeq = ''.join([ s for s in synBlockObj.quearySeq if s != \"-\" ])\n",
    "    if synBlockObj.orientation == \"-\":\n",
    "        qSeq = revCompSatsuma(qSeq)\n",
    "    qSeqLength,success,changed = len(qSeq),0,0\n",
    "    ########################\n",
    "    for i in range(0,howFar):\n",
    "        beginF = synBlockObj.quearyBegin+i\n",
    "        beginR = synBlockObj.quearyBegin-i\n",
    "        refSeqF = quearyRefDict[synBlockObj.quearyName][beginF:beginF+qSeqLength]\n",
    "        refSeqR = quearyRefDict[synBlockObj.quearyName][beginR:beginR+qSeqLength]\n",
    "        if qSeq == refSeqF:\n",
    "            success = 1\n",
    "            synBlockObj.quearyBegin = beginF\n",
    "            synBlockObj.quearyEnd = beginF+qSeqLength\n",
    "            if i != 0:\n",
    "                changed += 1\n",
    "            break\n",
    "        #####################\n",
    "        elif qSeq == refSeqR:\n",
    "            success = 1\n",
    "            synBlockObj.quearyBegin = beginR\n",
    "            synBlockObj.quearyEnd = beginR+qSeqLength\n",
    "            if i != 0:\n",
    "                changed += 1\n",
    "            break\n",
    "    ################\n",
    "    if success == 0:\n",
    "        print(\"### ERROR ###\")\n",
    "        print(synBlockObj.quearyName,synBlockObj.quearyBegin,synBlockObj.quearyEnd)\n",
    "        print(synBlockObj.length,synBlockObj.pID,synBlockObj.orientation)\n",
    "        return synBlockObj,changed\n",
    "    else:\n",
    "        return synBlockObj,changed\n",
    "    \n",
    "def readXHomology_toSynBlocks(XFile):\n",
    "    '''Takes in the MergeXCorrMatches.chained.out produced by Satsuma as input.\n",
    "       A SynBlock object is created for each sytenic pair found by satsuma and placed in\n",
    "       interval trees pertaining to the queary/target genomes. It's worth noting that the \n",
    "       target and queary interval trees both point to the same SynBlock object.\n",
    "       ############################## \n",
    "       XFile = filePath/To/MergeXFile\n",
    "       Output = [ synBlockObject , ... ] '''\n",
    "    inFile = open(XFile,'r')\n",
    "    ### Read header lines ###\n",
    "    for i in range(0,12):\n",
    "        inFile.readline()\n",
    "    ##############\n",
    "    synBlocks = []\n",
    "    synBlock = \"Dummy variable for the initialization of the alg... Is overwritten later. Not to worry\"\n",
    "    alignment = [ [],[],[] ]\n",
    "    aLCount = 1\n",
    "    while True:\n",
    "        lineMan = inFile.readline().strip('\\r').strip('\\n')\n",
    "        lineInfo = lineMan.split(' ')[0]\n",
    "        #######################\n",
    "        if lineInfo == \"Found\":\n",
    "            if len(alignment[0]) != 0:\n",
    "                synBlock.getSequences(alignment)\n",
    "                synBlocks.append(synBlock)\n",
    "            aLCount += 1\n",
    "            alignment = [ [],[],[] ]\n",
    "            infoLine = inFile.readline()\n",
    "            infoLine2 = inFile.readline()\n",
    "            inFile,synBlock = parseIntervalLines(inFile,infoLine,infoLine2 )   \n",
    "        #########################\n",
    "        elif lineInfo == \"Query\":\n",
    "            if len(alignment[0]) != 0:\n",
    "                synBlock.getSequences(alignment)\n",
    "                synBlocks.append(synBlock)\n",
    "            aLCount += 1\n",
    "            alignment = [ [],[],[] ]\n",
    "            infoLine2 = inFile.readline()\n",
    "            inFile,synBlock = parseIntervalLines(inFile,lineMan,infoLine2 )   \n",
    "        ###################\n",
    "        elif lineMan == '':\n",
    "            e1 = inFile.readline().strip('\\r').strip('\\n') ### Link line\n",
    "            e2 = inFile.readline().strip('\\r').strip('\\n') ### Sequence2 line\n",
    "            e3 = inFile.readline().strip('\\r').strip('\\n') ### Spacer line\n",
    "            e4 = inFile.readline().strip('\\r').strip('\\n') ### \"Protein coding line\"          \n",
    "            if \"Total\" not in e1:\n",
    "                alignment[1].append(e1)\n",
    "                alignment[2].append(e2)\n",
    "                inFile.readline()\n",
    "                inFile.readline()\n",
    "                \n",
    "            else:\n",
    "                print(\"End of file\")\n",
    "                print(\"Printing satsuma stats (last few lines of the MergeXHomology file)...\")\n",
    "                synBlock.getSequences(alignment)\n",
    "                synBlocks.append(synBlock)\n",
    "                aLCount += 1\n",
    "                alignment = [ [],[],[] ]\n",
    "                break\n",
    "        #####\n",
    "        else:\n",
    "            alignment[0].append(lineMan)\n",
    "            inFile,alignment = processAlignment(inFile,alignment)\n",
    "        ##########################\n",
    "        if aLCount % 1000000 == 0:\n",
    "            print(\"Alignments processed \"+str(aLCount))\n",
    "    ##############\n",
    "    inFile.close()\n",
    "    #########\n",
    "    print(e1)\n",
    "    print(e2)\n",
    "    print(e3)\n",
    "    print(e4)\n",
    "    print(\"Alignments put into interval tress \"+str(aLCount))\n",
    "    ################\n",
    "    return synBlocks\n",
    "\n",
    "def synBlocks_toIntervalTrees(synBlocks,qRef,tRef,howFar=2000):\n",
    "    ### howFar = Maximum distance to search in the quearyReference genome left/right of the quearyBegin position \n",
    "    altered = 0\n",
    "    qTrees,tTrees = {},{}\n",
    "    for i,sObj in enumerate(synBlocks):\n",
    "        if sObj.quearyName not in qTrees:\n",
    "            qTrees.update({ sObj.quearyName:intervaltree.IntervalTree() })\n",
    "        \n",
    "        if sObj.targetName not in tTrees:\n",
    "            tTrees.update({ sObj.targetName:intervaltree.IntervalTree() })\n",
    "        ##################################################################\n",
    "        \n",
    "        synBlock,c = preProcessIntervalCoords(sObj,qRef,tRef,howFar=howFar)\n",
    "        altered += c\n",
    "        \n",
    "        qTrees[sObj.quearyName].addi(sObj.quearyBegin,sObj.quearyEnd,sObj)\n",
    "        tTrees[sObj.targetName].addi(sObj.targetBegin,sObj.targetEnd,sObj)\n",
    "        \n",
    "        if i % 500000 == 0:\n",
    "            print(\"Synblocks procssed \"+str(i)+\"...\")\n",
    "    ###########################################################\n",
    "    print(\"Syntenic blocks with adjusted coords \"+str(altered))\n",
    "    return qTrees,tTrees,synBlocks       "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "**The general flow of the program is as follows**\n",
    "\n",
    "\n",
    "1. Read the queary and target reference genomes into memory in the form of { scaffoldName:scaffoldSequence } (one hashTable for each genome)  \n",
    "2. Read the satsuma MergeXHomology file into memory by storing the information in the form of a SynBlock Object. This is just an easy way to store/access all the information  \n",
    "3. \"Pre-process\" the SyntenyBlock Objects, by changing the quearyBegin and quearyEnd attributes if the satsuma coordinates are wrong.  \n",
    "We do this by matching up the quearySequence (output by satsuma) to the respective reference genome sequence, and pull the coordinates (if they have changed)  \n",
    "This is only an issue for the Queary genome that was used in the Satsuma run that generated all the data.  \n",
    "4. Add the queary and target information to their respective intervalTrees for a fast lookup  \n",
    "5. Any given queary/target genome position is now quearyable to obtain its sytenic partner"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-10-21T13:54:55.629135Z",
     "start_time": "2020-10-21T13:54:39.583183Z"
    },
    "hidden": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time to read fasta sequences to memory 16.03819990158081 seconds\n"
     ]
    }
   ],
   "source": [
    "startTime = time.time()\n",
    "marmDict = readFastaToMem(\"/fsimb/groups/imb-baumanngr/ao/HYBRID_COORDINATES/a_marmorata_renamed.fasta\")\n",
    "gulDict = readFastaToMem(\"/fsimb/groups/imb-baumanngr/ao/HYBRID_COORDINATES/a_gularis_unmasked_fixed.fasta\")\n",
    "print(\"Time to read fasta sequences to memory \"+str(time.time()-startTime)+\" seconds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-10-21T13:55:43.657519Z",
     "start_time": "2020-10-21T13:54:55.633449Z"
    },
    "hidden": true,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Alignments processed 1000000\n",
      "Alignments processed 1000000\n",
      "Alignments processed 2000000\n",
      "Alignments processed 2000000\n",
      "Alignments processed 2000000\n",
      "Alignments processed 2000000\n",
      "Alignments processed 2000000\n",
      "Alignments processed 2000000\n",
      "Alignments processed 2000000\n",
      "Alignments processed 2000000\n",
      "Alignments processed 2000000\n",
      "Alignments processed 2000000\n",
      "Alignments processed 2000000\n",
      "Alignments processed 2000000\n",
      "Alignments processed 2000000\n",
      "Alignments processed 2000000\n",
      "Alignments processed 2000000\n",
      "Alignments processed 2000000\n",
      "Alignments processed 2000000\n",
      "Alignments processed 2000000\n",
      "Alignments processed 2000000\n",
      "End of file\n",
      "Printing satsuma stats (last few lines of the MergeXHomology file)...\n",
      "Total bases aligned:               1230755977\n",
      "Total bases in target sequence:    1498801201\n",
      "Covered by at least one alignment: 1199035106 (79.9996 %)\n",
      "All bases in aligns: 2016206081\n",
      "Alignments put into interval tress 2311565\n",
      "Time to read satsuma alignments to syntenic blocks 48.01776909828186 seconds\n"
     ]
    }
   ],
   "source": [
    "startTime = time.time()\n",
    "synBlocks = readXHomology_toSynBlocks(\"/fsimb/groups/imb-baumanngr/ao/HYBRID_COORDINATES/MergeXCorrMatches.chained.out\")\n",
    "print(\"Time to read satsuma alignments to syntenic blocks \"+str(time.time()-startTime)+\" seconds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-10-21T14:00:33.432295Z",
     "start_time": "2020-10-21T13:55:43.659128Z"
    },
    "hidden": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Synblocks procssed 0...\n",
      "Synblocks procssed 500000...\n",
      "Synblocks procssed 1000000...\n",
      "Synblocks procssed 1500000...\n",
      "Synblocks procssed 2000000...\n",
      "Syntenic blocks with adjusted coords 59370\n",
      "Time to adjust coordinates and build interval tress 289.77008605003357 seconds\n"
     ]
    }
   ],
   "source": [
    "startTime = time.time()\n",
    "Q,T,synBlocks = synBlocks_toIntervalTrees(synBlocks,qRef=marmDict,tRef=gulDict,howFar=2000)\n",
    "print(\"Time to adjust coordinates and build interval tress \"+str(time.time()-startTime)+\" seconds\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Coordinate shift code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-09-22T16:30:52.155683Z",
     "start_time": "2020-09-22T16:30:52.153704Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": [
    "### From the satsuma manual/documentation \n",
    "### http://satsuma.sourceforge.net/faq.html\n",
    "\n",
    "### in principle, the TARGET should be the more COMPLETE genome\n",
    "\n",
    "### Satsuma chained out format = \n",
    "### Target sequence name (provided by fasta)\n",
    "### First target base\n",
    "### Last target base\n",
    "### Query sequence name (provided by fasta)\n",
    "### First query base\n",
    "### Last query base\n",
    "### Identity orientation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-10-21T14:01:05.141075Z",
     "start_time": "2020-10-21T14:01:05.115288Z"
    },
    "code_folding": [
     0,
     8,
     17,
     84
    ],
    "hidden": true
   },
   "outputs": [],
   "source": [
    "def printInfo(sObj):\n",
    "    print(sObj.quearyName,sObj.quearyBegin,sObj.quearyEnd)\n",
    "    print(sObj.targetName,sObj.targetBegin,sObj.targetEnd)\n",
    "    print(\"Orientation ( \"+str(sObj.orientation)+\" )\")\n",
    "    print(\"PercentIdentity:\"+str(sObj.pID))\n",
    "    print(\"Length:\"+str(sObj.length))\n",
    "    print(\"#####\")\n",
    "\n",
    "def printAln(sObj,step=100):\n",
    "    for i in range(0,sObj.length,step):\n",
    "        print(sObj.quearySeq[i:i+step])\n",
    "        print(sObj.linkSeq[i:i+step])\n",
    "        print(sObj.targetSeq[i:i+step])\n",
    "        print()\n",
    "        print(''.join([\"-\"]*step))\n",
    "        print()\n",
    "\n",
    "def quearySynBlock(synBlockObj,chrom,point,which=\"queary\"):\n",
    "    if synBlockObj.orientation == \"+\":\n",
    "        qCoord,tCoord = synBlockObj.quearyBegin,synBlockObj.targetBegin\n",
    "        successIndex = \"NA\"\n",
    "        for i in range(0,len(synBlockObj.quearySeq)):\n",
    "            if synBlockObj.quearySeq[i] != \"-\":\n",
    "                qCoord += 1\n",
    "                if which == \"queary\":\n",
    "                    if qCoord == point:\n",
    "                        successIndex = i\n",
    "            ###################################\n",
    "            if synBlockObj.targetSeq[i] != \"-\":\n",
    "                tCoord += 1\n",
    "                if which == \"target\":\n",
    "                    if tCoord == point:\n",
    "                        successIndex = i\n",
    "            ########################\n",
    "            if successIndex != \"NA\":\n",
    "                break\n",
    "        ########################\n",
    "        #if successIndex == \"NA\":\n",
    "        #    return \"NA\",\"NA\",\"NA\"\n",
    "        ##########################################################################################\n",
    "        ### Returns [ quearySeq,link,targetSeq ],[quearyName,quearyCoord],[targetName,targetCoord]\n",
    "        ###############################################\n",
    "        return [ synBlockObj.quearySeq[successIndex], \\\n",
    "                 synBlockObj.linkSeq[successIndex],\\\n",
    "                 synBlockObj.targetSeq[successIndex],\n",
    "                 synBlockObj.orientation ], \\\n",
    "               [synBlockObj.quearyName,qCoord],\\\n",
    "               [synBlockObj.targetName,tCoord]\n",
    "    ####################################\n",
    "    ####################################\n",
    "    ####################################\n",
    "    elif synBlockObj.orientation == \"-\":\n",
    "        if which == \"queary\":\n",
    "            difference = synBlockObj.quearyEnd - point + 1\n",
    "        elif which == \"target\":\n",
    "            difference = point - synBlockObj.targetBegin\n",
    "        #############################################################\n",
    "        qCoord,tCoord = synBlockObj.quearyEnd,synBlockObj.targetBegin\n",
    "        successIndex = \"NA\"\n",
    "        for i in range(0,len(synBlockObj.quearySeq)):\n",
    "            if synBlockObj.quearySeq[i] != \"-\":\n",
    "                qCoord -= 1\n",
    "                if which == \"queary\":\n",
    "                    difference -= 1\n",
    "            ###################################  \n",
    "            if synBlockObj.targetSeq[i] != \"-\":\n",
    "                tCoord += 1\n",
    "                if which == \"target\":\n",
    "                    difference -=1\n",
    "            ###################\n",
    "            if difference == 0:\n",
    "                successIndex = i\n",
    "                break\n",
    "        ########################\n",
    "        #if successIndex == \"NA\":\n",
    "        #    return \"NA\",\"NA\",\"NA\"   \n",
    "        ###############################################\n",
    "        return [ synBlockObj.quearySeq[successIndex], \\\n",
    "                 synBlockObj.linkSeq[successIndex],\\\n",
    "                 synBlockObj.targetSeq[successIndex],\n",
    "                 synBlockObj.orientation ], \\\n",
    "               [synBlockObj.quearyName,qCoord+1],\\\n",
    "               [synBlockObj.targetName,tCoord]\n",
    "\n",
    "def quearyTree(chrom,point,treeDict,which=\"queary\",method=\"first\"):\n",
    "    ### chrom,point = chromosome/position to look up in the interval tree\n",
    "    ### treeDict = { chromosomeName:intervalTree() }\n",
    "    ### which = \"queary\" or \"target\". This refers to the identity of the genome that the point belongs to. \n",
    "    ###         \"queary\" means the genome was used as the queary in the satsuma run. \n",
    "    ###         \"target\" means the genome was used as the target in the satsuma run. \n",
    "    \n",
    "    ### This will return a few different things...\n",
    "    ### syntenic coordinate found (or lack there of)\n",
    "    \n",
    "    ### The way Satsuma reports the information is with respect to the target genome that was used as input\n",
    "    ### Therefore, all coordinates and sequences for the target genome are with respect\n",
    "    ### to the (+) strand. The target genome acts like the reference genome for alignment\n",
    "    ### So for any given alignment that is reported, the target sequence is ALWAYS 5'-3' positive strand with repsect to the target reference genome\n",
    "    ### For an alignment where the orientation is (+), the queary sequence is also 5'-3' positive strand with repsect to the queary reference genome\n",
    "    ### For an alignment where the orientation is (-), the queary sequence is RT with respect to the queary reference genome.\n",
    "    ### The coordinates for both queary and target are again always with respect to the positive strand of the repsective reference genome.\n",
    "    ### It is only the queary sequence that will need to be reverseTranscribed with repsect to its reference genome (but the coords remain the same).\n",
    "    \n",
    "    \n",
    "    ### The information that is returned is in the form of \n",
    "    ### [ quearySequence,linkSequence,targetSequence,blockOrientation ] , [quearyName,quearyPoint] , [targetName,targetPoint]\n",
    "    ### We return both the queary and target information so that a single function can be used to investigate the sytenny\n",
    "    ### rather than two different ones for each queary/target genome\n",
    "    ### Will return \"NA\",\"NA\",\"NA\" if no sytenic coordinate is found in the other genome\n",
    "    \n",
    "    ### Generate a list of syntenic blocks that this point overlaps ###\n",
    "   \n",
    "    synBlockIntervals = list(treeDict[chrom][point-1])\n",
    "    blockCount = len(synBlockIntervals)\n",
    "    ###################\n",
    "    if blockCount == 1:\n",
    "        a,q,t = quearySynBlock(synBlockIntervals[0].data,chrom,point,which=which)\n",
    "    #####################\n",
    "    elif blockCount == 0:\n",
    "        a,q,t = \"NA\",\"NA\",\"NA\"\n",
    "    #######################\n",
    "    elif method == \"first\":\n",
    "        a,q,t = quearySynBlock(synBlockIntervals[0].data,chrom,point,which=which)\n",
    "    else:\n",
    "        ### If there are more than one syntenicBlocks overlapping the position then we return \"NA\"\n",
    "        ### This is the point in the algorithm in which a user defined function to determine the consensus sequence\n",
    "        ### could be utilized or some other method to deal with these edge cases. I have written a function to come\n",
    "        ### up with a consensus sequence, but how I define the consensus might be different from the users, so I\n",
    "        ### have left this part open to interpretation. ~1.5% of randomly quearied positions appear to have multiple\n",
    "        ### blocks overlapping them.\n",
    "        a,q,t = \"NA\",\"NA\",\"NA\"\n",
    "    ############\n",
    "    return a,q,t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "### Usage description/example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "The program is written for 1 based indices. Any position in IGV that you queary should come back exactly the same if you are using the quearyTree function.  \n",
    "Under the hood the quearyTree function will manually subtract 1 from the position, and use this to look up the overlapping sytnenicBlock(s). Then the original position will then be passed into the quearySynBlock function. \n",
    "\n",
    "But if the quearyTree function is not used, then any IGV coordinate you queary manually will come back as the position+1. So manually subtracting 1 from the original point is necessary when looking up the overlapping syntenicBlock(s), then the original point is used to obtain the coordinate.\n",
    "\n",
    "**Using the quearyTree function**  \n",
    "a,q,t = quearyTree(chrom,point,treeDict,which=\"queary\",method=\"first\")\n",
    "\n",
    "**Manually quearying the intervalTree**  \n",
    "blockMan = list(Q[chrom][**point-1**])[0].data  \n",
    "a,q,t = quearySynBlock(blockMan,chrom,**point**,which=\"queary\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-10-01T11:47:12.450935Z",
     "start_time": "2020-10-01T11:47:12.442214Z"
    },
    "hidden": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['T', ' ', 'G', '+'] ['Chromosome_1_marmoratus', 128638491] ['ScHqM4t_1_gularis', 26917]\n",
      "['T', ' ', 'G', '+'] ['Chromosome_1_marmoratus', 128638491] ['ScHqM4t_1_gularis', 26917]\n"
     ]
    }
   ],
   "source": [
    "blockMan = list(Q[\"Chromosome_1_marmoratus\"][128638491-1])[0].data\n",
    "a,q,t = quearySynBlock(blockMan,\"Chromosome_1_marmoratus\",128638491,which=\"queary\")\n",
    "print(a,q,t)\n",
    "a,q,t = quearyTree(\"Chromosome_1_marmoratus\",128638491,Q,which=\"queary\",method=\"first\")\n",
    "print(a,q,t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Recommended Usage with Pickle"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To avoid redoing the computational work needed to create the tree each time it's needed, it is recommended to save it as a pickled object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pickle.dump(query_tree, open(\"query_tree.p\", \"wb\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The tree can then be quickly loaded into memory and queried as needed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tree = pickle.load(open(\"query_tree.p\", \"rb\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "### Testing the code "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-09-22T13:37:41.233189Z",
     "start_time": "2020-09-22T13:37:41.227463Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": [
    "### Dumb/elegant way of assigning values to keys in a hashtable if the key doesn't exist ###\n",
    "### Essentially done in a single line as oppossed to my old method of querying the hashtable twice, resulting in more lines and slower runtime ###\n",
    "#counts = {}\n",
    "#for i in range(0,100):\n",
    "#    counts[1] = counts.get(1,0)+1\n",
    "#print(counts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "### For which==\"queary\" ###\n",
    "which = \"queary\"\n",
    "howMany = 100000\n",
    "fail,success = 0,0\n",
    "\n",
    "for c in marmDict:\n",
    "    if c not in Q:\n",
    "        continue\n",
    "    for i in range(1,howMany):\n",
    "        p = random.randint(1,len(marmDict[c]))\n",
    "        seq,qInfo,tInfo = quearyTree(c,p,Q,which=which,metho=\"first\")\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-09-22T15:01:26.848499Z",
     "start_time": "2020-09-22T14:56:39.190359Z"
    },
    "hidden": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['G', '|', 'G', '-'] 54 100.0\n",
      "['C', '|', 'C', '+'] 54 100.0\n",
      "['C', '|', 'C', '+'] 258 75.1938\n",
      "['G', '|', 'G', '-'] 213 66.1972\n",
      "['C', '|', 'C', '+'] 25 100.0\n",
      "['G', '|', 'G', '-'] 287 78.0488\n",
      "['G', '|', 'G', '-'] 54 100.0\n",
      "['G', '|', 'G', '-'] 157 53.5032\n",
      "###############\n",
      "###############\n",
      "###############\n",
      "###############\n",
      "###############\n",
      "['G', '|', 'G', '-'] 54 100.0\n",
      "['C', '|', 'C', '+'] 54 100.0\n",
      "['C', '|', 'C', '+'] 258 75.1938\n",
      "['G', '|', 'G', '-'] 213 66.1972\n",
      "['C', '|', 'C', '+'] 25 100.0\n",
      "['G', '|', 'G', '-'] 287 78.0488\n",
      "['G', '|', 'G', '-'] 54 100.0\n",
      "['G', '|', 'G', '-'] 157 53.5032\n",
      "###############\n",
      "###############\n",
      "###############\n",
      "###############\n",
      "###############\n",
      "['C', '|', 'C', '-'] 225 79.5556\n",
      "['C', ' ', 'G', '-'] 202 69.3069\n",
      "['C', ' ', '-', '-'] 182 64.8352\n",
      "['C', ' ', 'T', '-'] 188 63.8298\n",
      "['C', '|', 'C', '-'] 175 64.0\n",
      "['C', '|', 'C', '-'] 214 65.8879\n",
      "['C', '|', 'C', '-'] 250 56.0\n",
      "['C', '|', 'C', '-'] 188 65.9574\n",
      "###############\n",
      "###############\n",
      "###############\n",
      "###############\n",
      "###############\n",
      "['T', '|', 'T', '-'] 225 79.5556\n",
      "['T', ' ', 'G', '-'] 202 69.3069\n",
      "['T', ' ', '-', '-'] 182 64.8352\n",
      "['T', '|', 'T', '-'] 188 63.8298\n",
      "['T', '|', 'T', '-'] 175 64.0\n",
      "['T', '|', 'T', '-'] 214 65.8879\n",
      "['T', '|', 'T', '-'] 250 56.0\n",
      "['T', '|', 'T', '-'] 188 65.9574\n",
      "###############\n",
      "###############\n",
      "###############\n",
      "###############\n",
      "###############\n",
      "['A', '|', 'A', '-'] 225 79.5556\n",
      "['A', '|', 'A', '-'] 202 69.3069\n",
      "['A', ' ', '-', '-'] 182 64.8352\n",
      "['A', ' ', 'T', '-'] 188 63.8298\n",
      "['A', '|', 'A', '-'] 175 64.0\n",
      "['A', '|', 'A', '-'] 214 65.8879\n",
      "['A', '|', 'A', '-'] 250 56.0\n",
      "['A', '|', 'A', '-'] 188 65.9574\n",
      "###############\n",
      "###############\n",
      "###############\n",
      "###############\n",
      "###############\n",
      "['G', '|', 'G', '-'] 225 79.5556\n",
      "['G', '|', 'G', '-'] 202 69.3069\n",
      "['G', ' ', '-', '-'] 182 64.8352\n",
      "['G', ' ', 'A', '-'] 188 63.8298\n",
      "['G', '|', 'G', '-'] 175 64.0\n",
      "['G', '|', 'G', '-'] 214 65.8879\n",
      "['G', '|', 'G', '-'] 250 56.0\n",
      "['G', '|', 'G', '-'] 188 65.9574\n",
      "###############\n",
      "###############\n",
      "###############\n",
      "###############\n",
      "###############\n",
      "['C', '|', 'C', '-'] 225 79.5556\n",
      "['C', ' ', '-', '-'] 202 69.3069\n",
      "['C', '|', 'C', '-'] 182 64.8352\n",
      "['C', '|', 'C', '-'] 188 63.8298\n",
      "['C', '|', 'C', '-'] 175 64.0\n",
      "['C', '|', 'C', '-'] 214 65.8879\n",
      "['C', '|', 'C', '-'] 250 56.0\n",
      "['C', '|', 'C', '-'] 188 65.9574\n",
      "###############\n",
      "###############\n",
      "###############\n",
      "###############\n",
      "###############\n",
      "['C', '|', 'C', '+'] 138 54.3478\n",
      "['C', '|', 'C', '+'] 141 57.4468\n",
      "['C', '|', 'C', '+'] 184 59.2391\n",
      "['C', '|', 'C', '+'] 151 54.9669\n",
      "['C', '|', 'C', '+'] 171 63.7427\n",
      "['C', '|', 'C', '+'] 165 60.6061\n",
      "['C', '|', 'C', '+'] 158 57.5949\n",
      "['C', '|', 'C', '+'] 158 65.1899\n",
      "###############\n",
      "###############\n",
      "###############\n",
      "###############\n",
      "###############\n",
      "['C', '|', 'C', '+'] 138 54.3478\n",
      "['C', '|', 'C', '+'] 141 57.4468\n",
      "['C', '|', 'C', '+'] 184 59.2391\n",
      "['C', '|', 'C', '+'] 151 54.9669\n",
      "['C', '|', 'C', '+'] 171 63.7427\n",
      "['C', '|', 'C', '+'] 165 60.6061\n",
      "['C', '|', 'C', '+'] 158 57.5949\n",
      "['C', '|', 'C', '+'] 158 65.1899\n",
      "###############\n",
      "###############\n",
      "###############\n",
      "###############\n",
      "###############\n",
      "['A', '|', 'A', '+'] 138 54.3478\n",
      "['A', ' ', 'T', '+'] 141 57.4468\n",
      "['A', ' ', 'T', '+'] 184 59.2391\n",
      "['A', ' ', 'T', '+'] 151 54.9669\n",
      "['A', ' ', 'T', '+'] 171 63.7427\n",
      "['A', ' ', 'T', '+'] 165 60.6061\n",
      "['A', ' ', 'T', '+'] 158 57.5949\n",
      "['A', ' ', 'C', '+'] 158 65.1899\n",
      "###############\n",
      "###############\n",
      "###############\n",
      "###############\n",
      "###############\n",
      "['C', ' ', 'T', '+'] 138 54.3478\n",
      "['C', ' ', 'G', '+'] 141 57.4468\n",
      "['C', ' ', 'T', '+'] 184 59.2391\n",
      "['C', ' ', 'T', '+'] 151 54.9669\n",
      "['C', '|', 'C', '+'] 171 63.7427\n",
      "['C', '|', 'C', '+'] 165 60.6061\n",
      "['C', ' ', 'T', '+'] 158 57.5949\n",
      "['C', '|', 'C', '+'] 158 65.1899\n",
      "###############\n",
      "###############\n",
      "###############\n",
      "###############\n",
      "###############\n",
      "['T', ' ', 'A', '+'] 138 54.3478\n",
      "['T', ' ', 'C', '+'] 141 57.4468\n",
      "['T', '|', 'T', '+'] 184 59.2391\n",
      "['T', ' ', 'C', '+'] 151 54.9669\n",
      "['T', ' ', 'C', '+'] 171 63.7427\n",
      "['T', '|', 'T', '+'] 165 60.6061\n",
      "['T', '|', 'T', '+'] 158 57.5949\n",
      "['T', '|', 'T', '+'] 158 65.1899\n",
      "###############\n",
      "###############\n",
      "###############\n",
      "###############\n",
      "###############\n",
      "{0: 6973745, 1: 15618643, 2: 395144, 3: 9855, 5: 535, 4: 1802, 6: 218, 8: 12, 7: 22, 9: 1}\n"
     ]
    }
   ],
   "source": [
    "### For which==\"queary\" ###\n",
    "which = \"queary\"\n",
    "howMany = 1000000\n",
    "fail,success = 0,0\n",
    "\n",
    "countDict = {}\n",
    "\n",
    "for c in marmDict:\n",
    "    if c not in Q:\n",
    "        continue\n",
    "    for i in range(1,howMany):\n",
    "        p = random.randint(1,len(marmDict[c]))\n",
    "        synBlockIntervals = list(Q[c][p-1])\n",
    "        countDict[len(synBlockIntervals)] = countDict.get(len(synBlockIntervals),0)+1\n",
    "        #if len(synBlockIntervals) == 8:\n",
    "        #    for interval in synBlockIntervals:\n",
    "        #        sObj = interval.data\n",
    "        #        printInfo(sObj)\n",
    "        #        print()\n",
    "        #        printAln(sObj,step=200)\n",
    "        #        print(\"###############\")\n",
    "        #        print(\"###############\")\n",
    "        #        print(\"###############\")\n",
    "        #        print(\"###############\")\n",
    "        #        print(\"###############\")\n",
    "\n",
    "        \n",
    "        if len(synBlockIntervals) == 8:\n",
    "            for interval in synBlockIntervals:\n",
    "                sObj = interval.data\n",
    "                a,q,t = quearySynBlock(sObj,c,p,which=which)\n",
    "                print(a,sObj.length,sObj.pID)\n",
    "            print(\"###############\")\n",
    "            print(\"###############\")\n",
    "            print(\"###############\")\n",
    "            print(\"###############\")\n",
    "            print(\"###############\")\n",
    "                \n",
    "        \n",
    "print(countDict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-10-01T11:39:05.074550Z",
     "start_time": "2020-10-01T11:34:29.419398Z"
    },
    "code_folding": [],
    "hidden": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Success 1548317\n",
      "Fail 0\n"
     ]
    }
   ],
   "source": [
    "### For which==\"queary\" ###\n",
    "which = \"queary\"\n",
    "howMany = 100000\n",
    "fail,success = 0,0\n",
    "\n",
    "for c in marmDict:\n",
    "    if c not in Q:\n",
    "        continue\n",
    "    for i in range(1,howMany):\n",
    "        p = random.randint(1,len(marmDict[c]))\n",
    "        synBlockIntervals = list(Q[c][p-1])\n",
    "        \n",
    "        \n",
    "        if len(synBlockIntervals) == 1:\n",
    "            sObj = synBlockIntervals[0].data\n",
    "\n",
    "            a,q,t = (quearySynBlock(sObj,c,p,which=which))\n",
    "            if a == \"NA\":\n",
    "                continue\n",
    "                \n",
    "            ###########################\n",
    "            if sObj.orientation == \"-\":\n",
    "                quearyTest = revCompSatsuma(a[0])\n",
    "            else:\n",
    "                quearyTest = a[0]\n",
    "            targetTest = a[2]\n",
    "            #################\n",
    "            \n",
    "\n",
    "            ### Test for the \"queary\" information of the syntenic pair ###\n",
    "            if q[1] != p:\n",
    "                print(\"### Queary coordinate ERROR ###\")\n",
    "                print(c,p,q[1])\n",
    "                printInfo(sObj)\n",
    "                print()\n",
    "                printAln(sObj)\n",
    "            \n",
    "            if quearyTest != marmDict[c][p-1]:\n",
    "                print(\"### Queary coordinate/Sequence ERROR ###\")\n",
    "                print(c,p,q[1])\n",
    "                printInfo(sObj)\n",
    "                print()\n",
    "                printAln(sObj)\n",
    "                fail += 1\n",
    "            \n",
    "            ### Test for the \"target\" information of the syntenic pair ###\n",
    "            if targetTest == \"-\":\n",
    "                continue\n",
    "            if targetTest != gulDict[t[0]][t[1]-1]:\n",
    "                \n",
    "                print(\"### Target coordinate/Sequence ERROR ###\")\n",
    "                print(c,p,t[1])\n",
    "                print(t)\n",
    "                print(a)\n",
    "                print(gulDict[t[0]][t[1]-1])\n",
    "                printInfo(sObj)\n",
    "                print()\n",
    "                printAln(sObj)\n",
    "                fail += 1\n",
    "            \n",
    "            \n",
    "            \n",
    "            else:\n",
    "                success += 1\n",
    "\n",
    "\n",
    "print(\"Success \"+str(success))\n",
    "print(\"Fail \"+str(fail))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-10-01T11:43:15.072037Z",
     "start_time": "2020-10-01T11:40:36.878784Z"
    },
    "hidden": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Success 987976\n",
      "Fail 0\n"
     ]
    }
   ],
   "source": [
    "### For which==\"target\" ###\n",
    "which = \"target\"\n",
    "howMany = 10000\n",
    "fail,success = 0,0\n",
    "\n",
    "for c in gulDict:\n",
    "    if c not in T:\n",
    "        continue\n",
    "    for i in range(1,howMany):\n",
    "        p = random.randint(1,len(gulDict[c]))\n",
    "        synBlockIntervals = list(T[c][p-1])\n",
    "        \n",
    "        \n",
    "        if len(synBlockIntervals) == 1:\n",
    "            sObj = synBlockIntervals[0].data\n",
    "\n",
    "            a,q,t = (quearySynBlock(sObj,c,p,which=which))\n",
    "            if a == \"NA\":\n",
    "                continue\n",
    "                \n",
    "            ###########################\n",
    "            if sObj.orientation == \"-\":\n",
    "                quearyTest = revCompSatsuma(a[0])\n",
    "            else:\n",
    "                quearyTest = a[0]\n",
    "            targetTest = a[2]\n",
    "            #################\n",
    "            \n",
    "            ### Test for the \"queary\" information of the syntenic pair ###\n",
    "            if t[1] != p:\n",
    "                print(\"### Target coordinate ERROR ###\")\n",
    "                print(c,p,t[1])\n",
    "                printInfo(sObj)\n",
    "                print(a)\n",
    "                print()\n",
    "                printAln(sObj)\n",
    "            \n",
    "            if targetTest != gulDict[c][p-1]:\n",
    "                print(\"### Target coordinate/Sequence ERROR ###\")\n",
    "                print(c,p,t[1])\n",
    "                printInfo(sObj)\n",
    "                print(a)\n",
    "                print()\n",
    "                printAln(sObj)\n",
    "                fail += 1\n",
    "            \n",
    "            ### Test for the \"target\" information of the syntenic pair ###\n",
    "            if quearyTest == \"-\":\n",
    "                continue\n",
    "            if quearyTest != marmDict[q[0]][q[1]-1]:\n",
    "                \n",
    "                print(\"### Queary coordinate/Sequence ERROR ###\")\n",
    "                print(c,p,q[1])\n",
    "                print(t)\n",
    "                print(a)\n",
    "                print(gulDict[q[0]][q[1]-1])\n",
    "                printInfo(sObj)\n",
    "                print()\n",
    "                prin1tAln(sObj)\n",
    "                fail += 1\n",
    "            \n",
    "            \n",
    "            \n",
    "            else:\n",
    "                success += 1\n",
    "\n",
    "\n",
    "print(\"Success \"+str(success))\n",
    "print(\"Fail \"+str(fail))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}