OlowoDeduperPseudo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Deduper Assignment\n",
    "\n",
    "__Problem:__ Given SAM files of uniquely mapped reads, develop an algorithm to detect and remove PCR duplicates based on the bitwise flag, left most chromosome position, and the CIGAR string. Since SAM format only gives left most position, if a pa\n",
    "rticular read is stranded(bitwise flag == x), then I must write a function to extrapolate the starting position(right) of the read from POS and CIGAR string. The last step to confirm a duplicate would be to check UMIs to make sure they match. Or would checking UMI's first be more efficient?\n",
    "Start with single ended reads, and expand to paired end. Be sure your algorithm adjusts for soft clipping. Start with known UMIs then expand to randomers. If duplicate is found, which read do you keep? Jason keeps the first one he encounters, Pete chooses a random one. Up to writer. Should duplicate of choice be passed as a flag? \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_QNAME_ contains UMIs\n",
    "\n",
    "__Input File__\n",
    "\n",
    "NS500451:154:HWKTMBGXX:1:11101:24260:1121:CTGTTCAC      0       2       76814284        36      71M     *       0       0      TCCACCACAATCTTACCATCCTTCCTCCAGACCACATCGCGTTCTTTGTTCAACTCACAGCTCAAGTACAA    6AEEEEEEAEEAEEEEAAEEEEEEEEEAEEAEEAAEE<EEEEEEEEEAEEEEEEEAAEEAAAEAEEAEAE/ MD:Z:71 NH:i:1  HI:i:1  NM:i:0  SM:i:36 XQ:i:40 X2:i:0  XO:Z:UU\n",
    "\n",
    "NS500451:154:HWKTMBGXX:1:11101:18996:1145:CTGTTCAC      0       2       76814284       36      40M1I30M        *       0       0       GTCTCTTAGTTTATTATAAACCAGCTTCATAGGCCACAGAGGAAAAAGGACTATATACATACAGCCTTTTG    6AEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEEE MD:Z:53G16      NH:i:1  HI:i:1  NM:i:2  SM:i:36 XQ:i:40 X2:i:0     XO:Z:UU\n",
    "\n",
    "NS500451:154:HWKTMBGXX:1:11101:25533:1187:GTTCACCT      16       2       76743835        36      71M     *       0       0      CTTGGTAACTTTCAGAGAATTAGTCACAACTTCTGAAGCAACCACAGTCCATGCAAGTCGACTGGTTTCTC    6AEEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEEEEEAEEEEEEE<EEEEEEEEEEEEEEEEEEEE MD:Z:71 NH:i:1  HI:i:1  NM:i:0  SM:i:36 XQ:i:40 X2:i:0  XO:Z:UU\n",
    "\n",
    "NS500451:154:HWKTMBGXX:1:11101:20566:1080:TATAGNGC      16       2       76743825       36      61M10S     *       0       0      CCCCAAACAAGTGGTCCCCAAAAGAGACAGCCTCAAAATGGTCTAAGAAGCTGGCATAAAGGTCAGGAAAA    6AEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEAEEEEEEEEEEEEEEEEE MD:Z:71 NH:i:1  HI:i:1  NM:i:0  SM:i:36 XQ:i:40 X2:i:0  XO:Z:UU\n",
    "\n",
    "__Output File__\n",
    "\n",
    "NS500451:154:HWKTMBGXX:1:11101:24260:1121:CTGTTCAC      0       2       76814284        36      71M     *       0       0      TCCACCACAATCTTACCATCCTTCCTCCAGACCACATCGCGTTCTTTGTTCAACTCACAGCTCAAGTACAA    6AEEEEEEAEEAEEEEAAEEEEEEEEEAEEAEEAAEE<EEEEEEEEEAEEEEEEEAAEEAAAEAEEAEAE/ MD:Z:71 NH:i:1  HI:i:1  NM:i:0  SM:i:36 XQ:i:40 X2:i:0  XO:Z:UU\n",
    "\n",
    "NS500451:154:HWKTMBGXX:1:11101:25533:1187:GTTCACCT      16       2       76743835        36      71M     *       0       0      CTTGGTAACTTTCAGAGAATTAGTCACAACTTCTGAAGCAACCACAGTCCATGCAAGTCGACTGGTTTCTC    6AEEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEEEEEAEEEEEEE<EEEEEEEEEEEEEEEEEEEE MD:Z:71 NH:i:1  HI:i:1  NM:i:0  SM:i:36 XQ:i:40 X2:i:0  XO:Z:UU\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "unexpected EOF while parsing (<ipython-input-2-75f8f974490e>, line 31)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;36m  File \u001b[0;32m\"<ipython-input-2-75f8f974490e>\"\u001b[0;36m, line \u001b[0;32m31\u001b[0m\n\u001b[0;31m    \u001b[0m\n\u001b[0m    ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m unexpected EOF while parsing\n"
     ]
    }
   ],
   "source": [
    "# Convert SAM file to BAM file using samtools sort. Then sort by QNAME, which contains UMI's. This groups all identical UMIs\n",
    "# together. -n flag\n",
    "# Samtools sort Sam file in order to put reads of similar positions next to each other\n",
    "\n",
    "def openFile(file):\n",
    "    \"\"\"Use this function to open your read file\"\"\"\n",
    "    file = open(\"\")\n",
    "    \n",
    "def getInfo(singleRead):\n",
    "    \"\"\" Use this function to isolate bitwise flag, POS, and CIGAR string from everything else in SAM file\n",
    "    store read in seperate variable\"\"\"\n",
    "    return(UMI, bFlag, pos, cigar)\n",
    "def extrapolate(CIGAR, bitwiseFlag, POS):\n",
    "    \"\"\"Use this function to extrapolate the starting position of a stranded read  \"\"\"\n",
    "    if bitwiseFlag != 16:\n",
    "        #Check if this is a stranded read\n",
    "        print(\"extrapolate function called on non-stranded read, recheck loop\")\n",
    "        break\n",
    "    leftStart = POS\n",
    "    rightStart = 0\n",
    "    #Parse through CIGAR string and add to right/subtract starting position, depending on representation\n",
    "    return(rightStart)\n",
    "            \n",
    "#These nested dictionaries hold POS as keys, and values are an inner dictionary in which UMIs are keys\n",
    "# and reads are values \n",
    "UMIs = {}\n",
    "rUMIs = {}\n",
    "#Open file you are reading. Files are sorted by QNAME\n",
    "file = open(\"sAMFile.sam\", 'rt') #Dont forget to close\n",
    "outFile = open(\"sAMFile_deDuped.sam\", 'a+')\n",
    "\n",
    "# Open sorted file up. Read in all reads that have same Q name(UMI), making sure to initialize(reset) dictionary\n",
    "# every set of UMIs. Seperate out strands since the reads will be stored in two seperate dictionaries.\n",
    "# Store POS as key and read as value in each respective dictionary. POS stored as key because python dict\n",
    "# keys are unique so, if a duplicate is encountered, it will be skipped. If the read is stranded, use the \n",
    "# extrapolate function to adjust for things like softclipping & indels and obtain the true starting POS. \n",
    "# Before moving on to the next set of UMIs, write out the unique reads to out file since the dictionary \n",
    "# will be reset. This helps with the computational intensity. Wash, rinse, repeat for next set of UMIss\n",
    "\n",
    "\n",
    "while True:\n",
    "    sFlag = False\n",
    "    read = file.readlines(\"Read in first read\")\n",
    "    UMI, bit, POS, cigar = getINFO(read)\n",
    "    # Check if data is on opposite strand\n",
    "    if bit % 4 == 0:\n",
    "        sFlag = True\n",
    "    if sFlag = True:\n",
    "        POS = extrapolate(cigar, bit, POS)\n",
    "        rUMIs[POS] = {}\n",
    "        rUMIs[POS][UMI] = read\n",
    "    else:\n",
    "        UMIs[POS] = {}\n",
    "        UMIs[POS][UMI] = read\n",
    "    \n",
    "    \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "str"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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
}