BAM2csv_extraction_LTC1_v4.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2022.07.04 iPython notebook script to extract TAPE sequences from 10X sequencing reads\n",
    "# Junhong Choi\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Importing in tools\n",
    "\n",
    "import gzip\n",
    "import subprocess\n",
    "import os,sys,csv,re\n",
    "import itertools\n",
    "from collections import Counter\n",
    "from collections import OrderedDict\n",
    "from operator import itemgetter\n",
    "from datetime import datetime\n",
    "from optparse import OptionParser,OptionGroup\n",
    "\n",
    "# Usefule modules\n",
    "import numpy as np\n",
    "import edlib\n",
    "\n",
    "\n",
    "# Plot tools \n",
    "import pandas as pd\n",
    "\n",
    "# Additional function\n",
    "def reverse_complement(seq):\n",
    "    return complement(seq[::-1])\n",
    "\n",
    "def complement(seq):\n",
    "    complement = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A', 'N': 'N'} \n",
    "    bases = list(seq) \n",
    "    bases = [complement[base] for base in bases] \n",
    "    return ''.join(bases)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Building one-off white-list dictionary for the TargetBC:\n",
    "bases = ['A','C','T','G']\n",
    "\n",
    "TargetBCs_correct = ['']\n",
    "with open('TargetBC_WL19.csv') as f:\n",
    "    reader = csv.reader(f)\n",
    "    for row in reader:\n",
    "        TargetBCs_correct.append(row[0])\n",
    "\n",
    "\n",
    "TargetBC_WL = {}\n",
    "for BC in TargetBCs_correct:\n",
    "    for ii in range(0,len(BC)):\n",
    "        for base in bases:\n",
    "            if ii == 0:\n",
    "                mismatch = base + BC[1:]\n",
    "            elif ii < len(BC):\n",
    "                mismatch = BC[0:ii] + base + BC[ii+1:]\n",
    "            elif ii == len(BC):\n",
    "                mismatch = BC[0:ii] + base\n",
    "            try:\n",
    "                TargetBC_WL[mismatch] = BC\n",
    "            except KeyError:\n",
    "                continue\n",
    "\n",
    "                \n",
    "TAPE1_sequence = 'TGATGGTGAGCACG'\n",
    "TAPE1_WL = []\n",
    "for ii in range(0,len(TAPE1_sequence)):\n",
    "    for base in bases:\n",
    "        if ii == 0:\n",
    "            mismatch = base + TAPE1_sequence[1:]\n",
    "        elif ii < len(TAPE1_sequence):\n",
    "            mismatch = TAPE1_sequence[0:ii] + base + TAPE1_sequence[ii+1:]\n",
    "        elif ii == len(TAPE1_sequence):\n",
    "            mismatch = TAPE1_sequence[0:ii] + base\n",
    "        \n",
    "        TAPE1_WL.append(mismatch)\n",
    "\n",
    "        \n",
    "\n",
    "InsertBC_WL = ['']\n",
    "with open('InsertBC_WL19.csv') as f:\n",
    "    reader = csv.reader(f)\n",
    "    for row in reader:\n",
    "        InsertBC_WL.append(row[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "50729\n"
     ]
    }
   ],
   "source": [
    "c = 0\n",
    "InsertBC_filtered = 0\n",
    "cell_TAPE = {}\n",
    "InsertBC_table = {}\n",
    "TargetBC_table = {}\n",
    "\n",
    "# Opening the BAM extract file\n",
    "# Has [cBC][mBC(TargetBC-5xTAPE)][n_reads][n_UMI][n_reads_filtered][n_UMI_filtered]\n",
    "with open('CellByTape_10X_runAll_bamExtract_t2.csv') as f:\n",
    "    reader = csv.reader(f)\n",
    "    for row in reader:\n",
    "        if c == 0:\n",
    "            c+=1\n",
    "            continue\n",
    "        if c > 0:\n",
    "            if int(row[3]) > 1: # confirming threshold for UMI\n",
    "                TargetBC = row[1][35:43]\n",
    "                try:\n",
    "                    TargetBC_table[TargetBC] += 1\n",
    "                except KeyError:\n",
    "                    TargetBC_table[TargetBC] = 1\n",
    "                try:\n",
    "                    TargetBC = TargetBC_WL[TargetBC]\n",
    "                    TAPE = row[1][60:]\n",
    "                    for TAPE_var in TAPE1_WL:\n",
    "                        TAPE = TAPE.replace(TAPE_var, ',') # Usual TAPE-1\n",
    "                    TAPE = TAPE.replace('TGATGGTGAGCGAAGGAAA','') # Frequent ending\n",
    "                    TAPE = TAPE.replace('GGATGATG','GGA,TGATG') # Edit from chimeric ending \n",
    "                    TAPE = TAPE[:TAPE.rfind(',')] # Frequent ending\n",
    "                    #TAPE = TAPE[:TAPE.find('TGATGGTGAGCCAA')] # Frequent ending\n",
    "                    \n",
    "                    #InsertBC WhiteList\n",
    "                    InsertBC_WL_flag = 0\n",
    "                    insert_idx = [_.start() for _ in re.finditer(',', TAPE)] \n",
    "                    prev_idx = 0\n",
    "                    for idx in insert_idx:\n",
    "                        if TAPE[prev_idx:idx] not in InsertBC_WL:\n",
    "                            InsertBC_WL_flag += 1\n",
    "                            InsertBC_filtered += 1\n",
    "                        try:\n",
    "                            InsertBC_table[TAPE[prev_idx:idx]] += 1\n",
    "                        except KeyError:\n",
    "                            InsertBC_table[TAPE[prev_idx:idx]] = 1\n",
    "                        prev_idx = idx+1\n",
    "                    if TAPE[prev_idx:] not in InsertBC_WL: #Taking care of the last InsertBC\n",
    "                        InsertBC_WL_flag = 1\n",
    "                        InsertBC_filtered += 1\n",
    "                    try:\n",
    "                        InsertBC_table[TAPE[prev_idx:]] += 1\n",
    "                    except KeyError:\n",
    "                        InsertBC_table[TAPE[prev_idx:]] = 1\n",
    "            \n",
    "                    # If all InsertBC matched the WhiteList,\n",
    "                    if InsertBC_WL_flag == 0:\n",
    "                        nEntries = TAPE.count(',')\n",
    "                        TAPE = TAPE + ','*(4-nEntries)\n",
    "                        cell_target = row[0]+','+TargetBC\n",
    "                        TAPE_read = [TAPE, row[3]]\n",
    "                        #print(row[0]+','+TargetBC+'-'+TAPE+','+str(row[3]))\n",
    "                        c += 1\n",
    "                        try:\n",
    "                            past_read = cell_TAPE[cell_target][1]\n",
    "                            if int(past_read) <= int(row[3]):\n",
    "                                cell_TAPE[cell_target] = TAPE_read\n",
    "                                \n",
    "\n",
    "                        except KeyError:\n",
    "                            cell_TAPE[cell_target] = TAPE_read\n",
    "                except KeyError:\n",
    "                    continue\n",
    "\n",
    "print(InsertBC_filtered)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('CellByTargetBC5xTAPE_runAll_processed_MOI19.csv', 'w', newline='') as f0:\n",
    "    f0.write('Cell,TargetBC,Site1,Site2,Site3,Site4,Site5,nUMI\\n')\n",
    "    for key,value in cell_TAPE.items():\n",
    "        f0.write(key+','+value[0]+','+str(value[1])+'\\n')\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('LTC1_MOI19_TargetBC_table.csv', 'w', newline='') as f0:\n",
    "    f0.write('TargetBC,count\\n')\n",
    "    for key,value in TargetBC_table.items():\n",
    "        f0.write(key+','+str(value)+'\\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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}