EZCris_V1_py3.10_Jul03.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "cbd5798e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The position of interest is shown as bold: AGCTGAGATTGTGCCACTGCAGTCCAGCCTGGGCGACAGAGCAAGACTCCATCTCAAAACAACAACAACTAAAAACGAACAACAACAACAACAAAAAACTGAGGCCAGGTGCAGTGGCTCACACCTCTAATCCCAGCAATTTGGGAGGCTGAGGTGAGAGGATCACTTCAGCTCAGCAGTTCGAGACCAGCCCAGGCAACACAGGGAGACAGACCCTGTTATATTGCAGAGAGACCCCATCTCCACAAAATATAAAAATATTAGTCAGATGTGGTGGCATCCCTGCCGTCCCAGCTACTCAGGAGGCTGAGACAGGAGGATCGCTTGAGCCCAGGAGGTCAAGGCTGCAGTGAGCTGTGATCATGCCACTGCACTCCAGCCTGGGCCACAGAGCTAGACCCTGTCTCTAAAATTTTTAGAGACCTTATCTCTAAAAATAAATTAAATAAATAAACCGGGAGCACCTACTTTTTCTT\u001b[1mT\u001b[0mTTCTTTTACTTTTTTTTTTTTTTTTTGGAGACAGGGTCTCTATCACCCAGGCTGGAGTGCGGGGGCATGATCTTGGTTCACTGCAGCCTCGACCCCTCAGGCTCAGGCAGTTCTCCCACCTTAGCCTCCCCAGTAGCTGGGACTACAGGCACATGCCACCATGCCCGGCTAATTTTGTCTTTTTTTTTTTTTTTGGTAGAGACAGGGTCTCACCATGCTCCTCAAAACTCCTGGACTCGAGAGATCCTCCTGCCTCGGCCTCCCAAAATGCTGGGATTACGGTGTGAGCCGCTGTGCCCGGCTATTTTATTTTTAAATGAATAAAAGCTGGAGCACCCAACTTTTTTGTTGTTGTTTTTCTGAGACAGAGTTTTGCTCGTCACCCAGGCTGGAGTGCAGCGGCGCAATCTCGGCCCACTGCAACCTCTGCCTCCCGGGTTCAAGCGATTCTCCTGCCTCAGCCTCCTGAGTAG\n",
      "You need sgRNA for antisense sequence.\n",
      "Here is reverse complement:                CTACTCAGGAGGCTGAGGCAGGAGAATCGCTTGAACCCGGGAGGCAGAGGTTGCAGTGGGCCGAGATTGCGCCGCTGCACTCCAGCCTGGGTGACGAGCAAAACTCTGTCTCAGAAAAACAACAACAAAAAAGTTGGGTGCTCCAGCTTTTATTCATTTAAAAATAAAATAGCCGGGCACAGCGGCTCACACCGTAATCCCAGCATTTTGGGAGGCCGAGGCAGGAGGATCTCTCGAGTCCAGGAGTTTTGAGGAGCATGGTGAGACCCTGTCTCTACCAAAAAAAAAAAAAAAGACAAAATTAGCCGGGCATGGTGGCATGTGCCTGTAGTCCCAGCTACTGGGGAGGCTAAGGTGGGAGAACTGCCTGAGCCTGAGGGGTCGAGGCTGCAGTGAACCAAGATCATGCCCCCGCACTCCAGCCTGGGTGATAGAGACCCTGTCTCCAAAAAAAAAAAAAAAAAGTAAAAGAA\u001b[1mA\u001b[0mAAGAAAAAGTAGGTGCTCCCGGTTTATTTATTTAATTTATTTTTAGAGATAAGGTCTCTAAAAATTTTAGAGACAGGGTCTAGCTCTGTGGCCCAGGCTGGAGTGCAGTGGCATGATCACAGCTCACTGCAGCCTTGACCTCCTGGGCTCAAGCGATCCTCCTGTCTCAGCCTCCTGAGTAGCTGGGACGGCAGGGATGCCACCACATCTGACTAATATTTTTATATTTTGTGGAGATGGGGTCTCTCTGCAATATAACAGGGTCTGTCTCCCTGTGTTGCCTGGGCTGGTCTCGAACTGCTGAGCTGAAGTGATCCTCTCACCTCAGCCTCCCAAATTGCTGGGATTAGAGGTGTGAGCCACTGCACCTGGCCTCAGTTTTTTGTTGTTGTTGTTGTTCGTTTTTAGTTGTTGTTGTTTTGAGATGGAGTCTTGCTCTGTCGCCCAGGCTGGACTGCAGTGGCACAATCTCAGCT\n",
      "sgRNA Sequence            Cell     Enzyme  Efficacy(%)      CFD specificity score    GC(%)    Tm(°C)  Hairpin?\n",
      "------------------------  ------------------------------  -----------------------  -------  --------  ----------\n",
      "GTAAAAGAAAAAGAAAAAGT NGG  mESC     ABE-CP  50.5                             0.556       20     43.4   True\n",
      "                          U2OS     ABE-CP  49.1\n",
      "                          HEK293T  ABE-CP  52.4\n",
      "AAAGAAAAAGTAGGTGCTCC NGG                                                    1           40     51.24  True\n",
      "----------------------------------------------------------------------------------------\n",
      "Here are restriction enzymes(if any) to test on wild and mutated sequences. \n",
      "This is the wild type sequence:    TAAAAGAAAAAGAAAAA\n",
      "This is the mutant type sequence:  TAAAAGAAGAAGAAAAA\n",
      "These enzymes will potentially cut the wild sequence, but not mutant sequence.\n",
      " \n",
      "These enzymes will potentially cut the mutant sequence, but not wild sequence.\n",
      "MboII GAAGA\n",
      "Cco11366VI GAAGAA\n",
      " \n",
      "----------------------------------------------------------------------------------------\n",
      "\n",
      "\n",
      "Here are 10 primer pair combinations for ARMS PCR, suggested by Primer1(http://primer1.soton.ac.uk/primer1.html)\n",
      "\n",
      "******************************OUTPUT 1******************************\n",
      "Forward inner primer (T allele):                Melting temperature\n",
      "450 ATAAACCGGGAGCACCTACTTTTTCGTT 477                    68\n",
      "\n",
      "Reverse inner primer (C allele):\n",
      "506 TCCAAAAAAAAAAAAAAAAAGTAAAAGCAG 477                  63\n",
      "\n",
      "Forward outer primer (5' - 3'):\n",
      "255 AAAATATTAGTCAGATGTGGTGGCATCC 282                    65\n",
      "\n",
      "Reverse outer primer (5' - 3'):\n",
      "668 AAAAAAAAAAAAGACAAAATTAGCCGGG 641                    65\n",
      "\n",
      "Product size for T allele: 220\n",
      "Product size for C allele: 252\n",
      "Product size of two outer primers: 414\n",
      "\n",
      "******************************OUTPUT 2******************************\n",
      "Forward inner primer (T allele):                Melting temperature\n",
      "450 ATAAACCGGGAGCACCTACTTTTTCGTT 477                    68\n",
      "\n",
      "Reverse inner primer (C allele):\n",
      "506 TCCAAAAAAAAAAAAAAAAAGTAAAAGCAG 477                  63\n",
      "\n",
      "Forward outer primer (5' - 3'):\n",
      "216 CTGTTATATTGCAGAGAGACCCCATCTC 243                    65\n",
      "\n",
      "Reverse outer primer (5' - 3'):\n",
      "668 AAAAAAAAAAAAGACAAAATTAGCCGGG 641                    65\n",
      "\n",
      "Product size for T allele: 220\n",
      "Product size for C allele: 291\n",
      "Product size of two outer primers: 453\n",
      "\n",
      "******************************OUTPUT 3******************************\n",
      "Forward inner primer (T allele):                Melting temperature\n",
      "450 ATAAACCGGGAGCACCTACTTTTTCGTT 477                    68\n",
      "\n",
      "Reverse inner primer (C allele):\n",
      "506 TCCAAAAAAAAAAAAAAAAAGTAAAAGCAG 477                  63\n",
      "\n",
      "Forward outer primer (5' - 3'):\n",
      "254 AAAAATATTAGTCAGATGTGGTGGCATCC 282                   65\n",
      "\n",
      "Reverse outer primer (5' - 3'):\n",
      "668 AAAAAAAAAAAAGACAAAATTAGCCGGG 641                    65\n",
      "\n",
      "Product size for T allele: 220\n",
      "Product size for C allele: 253\n",
      "Product size of two outer primers: 415\n",
      "\n",
      "******************************OUTPUT 4******************************\n",
      "Forward inner primer (T allele):                Melting temperature\n",
      "450 ATAAACCGGGAGCACCTACTTTTTCGTT 477                    68\n",
      "\n",
      "Reverse inner primer (C allele):\n",
      "506 TCCAAAAAAAAAAAAAAAAAGTAAAAGCAG 477                  63\n",
      "\n",
      "Forward outer primer (5' - 3'):\n",
      "258 ATATTAGTCAGATGTGGTGGCATCCCT 284                     65\n",
      "\n",
      "Reverse outer primer (5' - 3'):\n",
      "668 AAAAAAAAAAAAGACAAAATTAGCCGGG 641                    65\n",
      "\n",
      "Product size for T allele: 220\n",
      "Product size for C allele: 249\n",
      "Product size of two outer primers: 411\n",
      "\n",
      "******************************OUTPUT 5******************************\n",
      "Forward inner primer (T allele):                Melting temperature\n",
      "450 ATAAACCGGGAGCACCTACTTTTTCGTT 477                    68\n",
      "\n",
      "Reverse inner primer (C allele):\n",
      "506 TCCAAAAAAAAAAAAAAAAAGTAAAAGCAG 477                  63\n",
      "\n",
      "Forward outer primer (5' - 3'):\n",
      "227 CAGAGAGACCCCATCTCCACAAAATAT 253                     65\n",
      "\n",
      "Reverse outer primer (5' - 3'):\n",
      "668 AAAAAAAAAAAAGACAAAATTAGCCGGG 641                    65\n",
      "\n",
      "Product size for T allele: 220\n",
      "Product size for C allele: 280\n",
      "Product size of two outer primers: 442\n",
      "\n",
      "******************************OUTPUT 6******************************\n",
      "Forward inner primer (T allele):                Melting temperature\n",
      "450 ATAAACCGGGAGCACCTACTTTTTCGTT 477                    68\n",
      "\n",
      "Reverse inner primer (C allele):\n",
      "506 TCCAAAAAAAAAAAAAAAAAGTAAAAGCAG 477                  63\n",
      "\n",
      "Forward outer primer (5' - 3'):\n",
      "218 GTTATATTGCAGAGAGACCCCATCTCC 244                     65\n",
      "\n",
      "Reverse outer primer (5' - 3'):\n",
      "668 AAAAAAAAAAAAGACAAAATTAGCCGGG 641                    65\n",
      "\n",
      "Product size for T allele: 220\n",
      "Product size for C allele: 289\n",
      "Product size of two outer primers: 451\n",
      "\n",
      "******************************OUTPUT 7******************************\n",
      "Forward inner primer (T allele):                Melting temperature\n",
      "450 ATAAACCGGGAGCACCTACTTTTTCGTT 477                    68\n",
      "\n",
      "Reverse inner primer (C allele):\n",
      "506 TCCAAAAAAAAAAAAAAAAAGTAAAAGCAG 477                  63\n",
      "\n",
      "Forward outer primer (5' - 3'):\n",
      "215 CCTGTTATATTGCAGAGAGACCCCATC 241                     65\n",
      "\n",
      "Reverse outer primer (5' - 3'):\n",
      "668 AAAAAAAAAAAAGACAAAATTAGCCGGG 641                    65\n",
      "\n",
      "Product size for T allele: 220\n",
      "Product size for C allele: 292\n",
      "Product size of two outer primers: 454\n",
      "\n",
      "******************************OUTPUT 8******************************\n",
      "Forward inner primer (T allele):                Melting temperature\n",
      "450 ATAAACCGGGAGCACCTACTTTTTCGTT 477                    68\n",
      "\n",
      "Reverse inner primer (C allele):\n",
      "506 TCCAAAAAAAAAAAAAAAAAGTAAAAGCAG 477                  63\n",
      "\n",
      "Forward outer primer (5' - 3'):\n",
      "253 TAAAAATATTAGTCAGATGTGGTGGCATCC 282                  65\n",
      "\n",
      "Reverse outer primer (5' - 3'):\n",
      "668 AAAAAAAAAAAAGACAAAATTAGCCGGG 641                    65\n",
      "\n",
      "Product size for T allele: 220\n",
      "Product size for C allele: 254\n",
      "Product size of two outer primers: 416\n",
      "\n",
      "******************************OUTPUT 9******************************\n",
      "Forward inner primer (T allele):                Melting temperature\n",
      "450 ATAAACCGGGAGCACCTACTTTTTCGTT 477                    68\n",
      "\n",
      "Reverse inner primer (C allele):\n",
      "506 TCCAAAAAAAAAAAAAAAAAGTAAAAGCAG 477                  63\n",
      "\n",
      "Forward outer primer (5' - 3'):\n",
      "259 TATTAGTCAGATGTGGTGGCATCCCT 284                      65\n",
      "\n",
      "Reverse outer primer (5' - 3'):\n",
      "668 AAAAAAAAAAAAGACAAAATTAGCCGGG 641                    65\n",
      "\n",
      "Product size for T allele: 220\n",
      "Product size for C allele: 248\n",
      "Product size of two outer primers: 410\n",
      "\n",
      "******************************OUTPUT 10******************************\n",
      "Forward inner primer (T allele):                Melting temperature\n",
      "450 ATAAACCGGGAGCACCTACTTTTTCGTT 477                    68\n",
      "\n",
      "Reverse inner primer (C allele):\n",
      "506 TCCAAAAAAAAAAAAAAAAAGTAAAAGCAG 477                  63\n",
      "\n",
      "Forward outer primer (5' - 3'):\n",
      "258 ATATTAGTCAGATGTGGTGGCATCCC 283                      65\n",
      "\n",
      "Reverse outer primer (5' - 3'):\n",
      "668 AAAAAAAAAAAAGACAAAATTAGCCGGG 641                    65\n",
      "\n",
      "Product size for T allele: 220\n",
      "Product size for C allele: 249\n",
      "Product size of two outer primers: 411\n",
      "\n",
      "\n",
      "Completed\n",
      "----------------------------------------------------------------------------------------\n",
      "\n",
      "\n",
      "******************************************************Disclaimer:***********************************************\n",
      "This code is for educational purposes only. Use it at your own risk. The author of this code makes no warranties\n",
      "about the accuracy or suitability of this code for any purpose. By using this code, you agree to hold the author\n",
      "harmless from any damages that may arise from its use.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "##V5: Jul 02, 13:00,  check for RFLP RE and CFD score \n",
    "\n",
    "\n",
    "\n",
    "########libraries ######################\n",
    "import re\n",
    "from Bio.Seq import Seq\n",
    "from tabulate import tabulate\n",
    "import primer3\n",
    "from Bio.Restriction import AllEnzymes\n",
    "import subprocess\n",
    "import os\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "######Functions #######################\n",
    "def validate_sequence(sequence):\n",
    "    \"\"\"\n",
    "    Validates the input DNA sequence.\n",
    "    \"\"\"\n",
    "    valid_chars = set(\"ACTG\")\n",
    "    if not set(sequence.upper()).issubset(valid_chars):\n",
    "        raise ValueError(\"Invalid DNA sequence. Only A, C, T, G characters are allowed.\")\n",
    "\n",
    "def validate_nucleotide_position(position, sequence):\n",
    "    \"\"\"\n",
    "    Validates the nucleotide position.\n",
    "    \"\"\"\n",
    "    if position < 1 or position > len(sequence):\n",
    "        raise ValueError(\"Invalid nucleotide position. Please enter a valid position within the sequence.\")\n",
    "\n",
    "def get_search_strand(sequence, position):\n",
    "    \"\"\"\n",
    "    Determines the search strand based on the nucleotide at the position of interest.\n",
    "    \"\"\"\n",
    "    nucleotide = sequence[position - 1]\n",
    "    return \"sense\" if nucleotide in [\"A\", \"C\"] else \"antisense\"\n",
    "\n",
    "def calculate_cg_percent(sequence):\n",
    "    \"\"\"\n",
    "    Calculates the CG percent of a given sequence.\n",
    "    \"\"\"\n",
    "    cg_count = sequence.count(\"C\") + sequence.count(\"G\")\n",
    "    return (cg_count / len(sequence)) * 100\n",
    "\n",
    "def find_pam_positions(sequence, search_strand):\n",
    "    \"\"\"\n",
    "    Finds all PAM positions on the selected strand.\n",
    "    \"\"\"\n",
    "    if search_strand == \"sense\":\n",
    "        pam_motif = re.compile(\"(?=AGG|TGG|CGG|GGG)\")\n",
    "    else:\n",
    "        pam_motif = re.compile(\"(?=AGG|TGG|CGG|GGG)\") #?=CCA|CCT|CCC|CCG #\n",
    "    return [match.start() for match in pam_motif.finditer(sequence)]\n",
    "\n",
    "\n",
    "\n",
    "def find_sorted_efficacies(BE_profile, nucleotide_pos):\n",
    "    efficacies = []\n",
    "\n",
    "    for cell, enzyme, *efficacy_values in BE_profile:\n",
    "        efficacy = efficacy_values[pos_grna]\n",
    "        nucleotide = sequence[nucleotide_pos - 1]\n",
    "        if nucleotide == 'A':\n",
    "            if enzyme in ['ABE', 'ABE-CP']:\n",
    "                if efficacy > 0:\n",
    "                    efficacies.append((cell, enzyme, efficacy))\n",
    "        else:\n",
    "            if enzyme in ['BE4', 'BE4-CP', 'evoA-BE4', 'AID-BE4', 'eA3A-BE4', 'CDA-BE4']:\n",
    "                if efficacy > 0:\n",
    "                    efficacies.append((cell, enzyme, efficacy))\n",
    "\n",
    "    sorted_efficacies = sorted(efficacies, reverse=True)\n",
    "    return sorted_efficacies\n",
    "\n",
    "\n",
    "def get_re_wild(wild, mut):\n",
    "    enzymes = AllEnzymes\n",
    "    for enzyme in enzymes:\n",
    "        if enzyme.site in wild and enzyme.site not in mut:\n",
    "            print(enzyme, enzyme.site)\n",
    "    print(\" \")\n",
    "\n",
    "def get_re_mut(wild, mut):\n",
    "    enzymes = AllEnzymes\n",
    "    for enzyme in enzymes:\n",
    "        if enzyme.site in mut and enzyme.site not in wild:\n",
    "            print(enzyme, enzyme.site)\n",
    "    print(\" \")\n",
    "\n",
    "def generate_mutant_sequence(sequence, nucleotide_pos):\n",
    "    pos9 = nucleotide_pos - 9\n",
    "    pos8 = nucleotide_pos + 8\n",
    "\n",
    "    wild = sequence[pos9:pos8]\n",
    "    #mut = wild[:8] + str(Seq(wild[8]).complement()) + wild[9:]\n",
    "    # if search_strand== \"sense\": \n",
    "    #     sequence=sequence\n",
    "    # else: \n",
    "    #     sequence= revcomp_sequence\n",
    "\n",
    "    #wild = sequence[pos9:pos8]\n",
    "    nuc_mut = Seq(wild[8])\n",
    "\n",
    "    if nuc_mut == \"A\":\n",
    "        mut = wild[:8] + \"G\" + wild[9:]\n",
    "    elif nuc_mut == \"C\":\n",
    "        mut = wild[:8] + \"T\" + wild[9:]\n",
    "    else:\n",
    "        mut = wild  # No mutation needed for other nucleotides\n",
    "\n",
    "    print(\"This is the wild type sequence:   \", wild)\n",
    "    print(\"This is the mutant type sequence: \", mut)\n",
    "\n",
    "    print(\"These enzymes will potentially cut the wild sequence, but not mutant sequence.\")\n",
    "    get_re_wild(wild, mut)\n",
    "    \n",
    "    \n",
    "    print(\"These enzymes will potentially cut the mutant sequence, but not wild sequence.\")\n",
    "    get_re_mut(wild, mut)\n",
    "\n",
    "####function for CFD score calculation \n",
    "def cfd_score_calculator(wt, off):\n",
    "    command = f\"python cfd_score_calculator.py --wt {wt} --off {off}\"\n",
    "    process = subprocess.Popen(command, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)\n",
    "    stdout, stderr = process.communicate()\n",
    "\n",
    "    # Check if the subprocess was executed successfully\n",
    "    if process.returncode != 0:\n",
    "        raise RuntimeError(f\"Error executing cfd_score_calculator.py: {stderr.decode()}\")\n",
    "\n",
    "    # Extract the numerical value from the output\n",
    "    output = stdout.decode().strip()\n",
    "    match = re.search(r\"(\\d+\\.\\d+)\", output)\n",
    "    if match:\n",
    "        score = float(match.group(1))\n",
    "    else:\n",
    "        raise ValueError(\"Could not extract CFD score from the output.\")\n",
    "\n",
    "    return round(score, 3)\n",
    "\n",
    "####\n",
    "def run_doPrimer(org_sequence, org_nucleotide_pos, org_nucleotide, mut_nucleotide):\n",
    "    # Change permissions of the shell script\n",
    "    os.chmod(\"run_doPrimer.sh\", 0o755)\n",
    "    \n",
    "    # Construct and execute the command\n",
    "    command = f\"./run_doPrimer.sh '{org_sequence} {org_nucleotide_pos} {org_nucleotide} {mut_nucleotide} 28 30 26 200 300 50 1.5 1.1 65 80 50 80 10 8.0 3.0 50 50 10'\"\n",
    "    subprocess.run(command, shell=True)\n",
    "\n",
    "\n",
    "######################################################################################################\n",
    "##################Base editor activity profiles#######################################################\n",
    "##################Reference: Mandana Arbab et al. Cell. 2020. doi:10.1016/j.cell.2020.05.037##########\n",
    "######################################################################################################\n",
    "BE_profile = [[\n",
    "  'mESC', 'ABE', 0.0, 0.0, 0.0, 37.4, 85.1, 100.0, 77.7, 43.1, 0.0, 0.0, 0.0,\n",
    "  0.0, 0.0, 0.0, 0.0, 0, 0, 0, 0, 0], [\n",
    "  'mESC', 'ABE-CP', 0.0, 0.0, 41.0, 83.9, 91.9, 100.0, 78.8, 63.1, 52.6, 50.5,\n",
    "  0.0, 0.0, 0.0, 0.0, 0.0, 0, 0, 0, 0, 0], [\n",
    "  'mESC', 'BE4', 0.0, 0.0, 31.8, 63.2, 90.3, 100.0, 87.0, 62.0, 31.4, 0.0, 0.0,\n",
    "  0.0, 0.0, 0.0, 0.0, 0, 0, 0, 0, 0], [\n",
    "  'mESC', 'BE4-CP', 0.0, 42.5, 47.9, 69.3, 80.8, 92.1, 97.0, 100.0, 97.3, 91.7,\n",
    "  76.8, 65.9, 59.6, 38.5, 34.7, 0, 0, 0, 0, 0], [\n",
    "  'mESC', 'evoA-BE4', 0.0, 44.4, 72.1, 97.3, 100.0, 98.9, 86.7, 63.6, 36.5, 0.0,\n",
    "  0.0, 0.0, 0.0, 0.0, 0.0, 0, 0, 0, 0, 0], [\n",
    "  'mESC', 'AID-BE4', 46.1, 62.4, 82.4, 87.2, 90.2, 95.9, 100.0, 86.6, 76.9, 56.5,\n",
    "  32.3, 0.0, 0.0, 0.0, 0.0, 0, 0, 0, 0, 0], [\n",
    "  'mESC', 'eA3A-BE4', 0.0, 0.0, 0.0, 41.6, 83.6, 100.0, 98.5, 75.7, 35.9, 0.0, 0.0,\n",
    "  0.0, 0.0, 0.0, 0.0, 0, 0, 0, 0, 0], [\n",
    "  'mESC', 'CDA-BE4', 88.4, 98.9, 97.8, 99.4, 100.0, 93.5, 90.4, 68.8, 40.1, 0.0,\n",
    "  0.0, 0.0, 0.0, 0.0, 0.0, 0, 0, 0, 0, 0], [\n",
    "  'HEK293T', 'ABE', 0.0, 0.0, 0.0, 54.2, 74.2, 100.0, 86.1, 50.2, 0.0, 0.0, 0.0, 0.0,\n",
    "  0.0, 0.0, 0.0, 0, 0, 0, 0, 0], [\n",
    "  'HEK293T', 'ABE-CP', 0.0, 0.0, 39.4, 59.4, 75.7, 77.4, 100.0, 76.1, 58.3, 52.4, 44.6, \n",
    "  0.0, 0.0, 0.0, 0.0, 0, 0, 0, 0, 0], [\n",
    "  'HEK293T', 'BE4', 0.0, 0.0, 36.1, 68.4, 100.0, 90.7, 72.2, 59.0, 37.3, 0.0, 0.0, 0.0,\n",
    "  0.0, 0.0, 0.0, 0, 0, 0, 0, 0], [\n",
    "  'HEK293T', 'BE4-CP', 0.0, 35.8, 50.4, 75.9, 81.3, 86.8, 83.7, 84.6, 100.0, 92.9, 92.0,\n",
    "  59.5, 56.1, 56.2, 37.9, 0, 0, 0, 0, 0], [\n",
    "  'HEK293T', 'evoA-BE4', 0.0, 39.7, 70.8, 97.7, 100.0, 94.8, 79.6, 67.1, 39.0, 0.0, 0.0,\n",
    "  0.0, 0.0, 0.0, 0.0, 0, 0, 0, 0, 0], [\n",
    "  'HEK293T', 'AID-BE4', 49.0, 62.0, 93.1, 87.9, 82.9, 100.0, 82.9, 97.0, 75.1, 65.6, 33.8,\n",
    "  33.6, 0.0, 0.0, 0.0, 0, 0, 0, 0, 0], [\n",
    "  'HEK293T', 'eA3A-BE4', 0.0, 0.0, 0.0, 35.8, 85.2, 100.0, 76.4, 73.7, 35.7, 0.0, 0.0, 0.0,\n",
    "  0.0, 0.0, 0.0, 0, 0, 0, 0, 0], [\n",
    "  'HEK293T', 'CDA-BE4', 80.8, 79.3, 84.6, 100.0, 83.3, 89.0, 80.0, 73.6, 57.3, 37.7, 0.0,\n",
    "  0.0, 0.0, 0.0, 0.0, 0, 0, 0, 0, 0], [\n",
    "  'HEK293T', 'eA3A-BE5', 0.0, 0.0, 0.0, 42.1, 67.9, 95.9, 100.0, 65.6, 32.1, 0.0, 0.0, 0.0,\n",
    "  0.0, 0.0, 0.0, 0, 0, 0, 0, 0], [\n",
    "  'U2OS', 'ABE', 0.0, 0.0, 0.0, 38.9, 100.0, 97.1, 83.5, 43.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\n",
    "  0.0, 0, 0, 0, 0, 0], [\n",
    "  'U2OS', 'ABE-CP', 0.0, 0.0, 34.1, 71.5, 83.0, 100.0, 90.1, 82.2, 55.1, 49.1, 44.4, 0.0, 0.0,\n",
    "  0.0, 0.0, 0, 0, 0, 0, 0], [\n",
    "  'U2OS', 'BE4', 0.0, 0.0, 41.1, 66.3, 100.0, 99.0, 93.2, 56.1, 40.3, 0.0, 0.0, 0.0, 0.0, 0.0,\n",
    "  0.0, 0, 0, 0, 0, 0], [\n",
    "  'U2OS', 'BE4-CP', 0.0, 34.5, 57.0, 64.9, 90.5, 82.7, 82.9, 100.0, 94.1, 94.9, 79.9, 81.0,\n",
    "  50.0, 32.0, 0.0, 0, 0, 0, 0, 0], [\n",
    "  'U2OS', 'evoA-BE4', 0.0, 32.1, 64.3, 92.4, 91.7, 100.0, 79.9, 61.2, 45.7, 0.0, 0.0, 0.0, 0.0,\n",
    "  0.0, 0.0, 0, 0, 0, 0, 0], [\n",
    "  'U2OS', 'AID-BE4', 49.0, 62.0, 93.1, 87.9, 82.9, 100.0, 82.9, 97.0, 75.1, 65.6, 33.8, 33.6,\n",
    "  0.0, 0.0, 0.0, 0, 0, 0, 0, 0], [\n",
    "  'U2OS', 'eA3A-BE4', 0.0, 0.0, 0.0, 0.0, 74.7, 100.0, 99.9, 65.7, 38.0, 0.0, 0.0, 0.0, 0.0, 0.0,\n",
    "  0.0, 0, 0, 0, 0, 0], [\n",
    "  'U2OS', 'CDA-BE4', 80.8, 79.3, 84.6, 100.0, 83.3, 89.0, 80.0, 73.6, 57.3, 37.7, 0.0, 0.0, 0.0,\n",
    "  0.0, 0.0, 0, 0, 0, 0, 0]]\n",
    "########################################################################################################\n",
    "\n",
    "\n",
    "# Get input sequence from user\n",
    "sequence = input(\"Please enter a DNA sequence(with at least 20nt flanking in upstream and downstream of position of the nucleotide of interest): \")\n",
    "validate_sequence(sequence)\n",
    "sequence = sequence.upper()\n",
    "org_sequence=sequence \n",
    "\n",
    "# Get nucleotide position of interest from user\n",
    "nucleotide_pos = int(input(\"Please enter the position of the nucleotide of interest: \"))\n",
    "validate_nucleotide_position(nucleotide_pos, sequence)\n",
    "\n",
    "\n",
    "org_nucleotide_pos=nucleotide_pos\n",
    "\n",
    "###name of nucleotide of interest \n",
    "org_nucleotide = sequence[org_nucleotide_pos - 1]\n",
    "\n",
    "###mutation type for ARMS PCR\n",
    "# Mutation type for ARMS PCR\n",
    "if org_nucleotide == \"A\":\n",
    "    mut_nucleotide = \"G\"\n",
    "elif org_nucleotide == \"C\":\n",
    "    mut_nucleotide = \"T\"\n",
    "elif org_nucleotide == \"T\":\n",
    "    mut_nucleotide = \"C\"\n",
    "elif org_nucleotide == \"G\":\n",
    "    mut_nucleotide = \"A\"\n",
    "else:\n",
    "    \n",
    "    mut_nucleotide = None  # This is a wrong nt letter. \n",
    "\n",
    "\n",
    "\n",
    "# Create a pattern to match the nucleotide of interest\n",
    "nucleotide_pattern = re.compile(f\"(?P<before>.{{0,{nucleotide_pos-1}}})(?P<nucleotide>.)(?P<after>.*)\")\n",
    "\n",
    "# Create a formatted string with the bold nucleotide position\n",
    "formatted_sequence = nucleotide_pattern.sub(\"\\g<before>\\033[1m\\g<nucleotide>\\033[0m\\g<after>\", sequence)\n",
    "\n",
    "# Print the formatted sequence\n",
    "print(f\"The position of interest is shown as bold: {formatted_sequence}\")\n",
    "\n",
    "# Determine the search strand based on the nucleotide at the position of interest\n",
    "search_strand = get_search_strand(sequence, nucleotide_pos)\n",
    "print(f\"You need sgRNA for {search_strand} sequence.\")\n",
    "\n",
    "\n",
    "\n",
    "# Reverse complement the sequence if it's on the antisense strand\n",
    "if search_strand == \"antisense\":\n",
    "    revcomp_sequence = str(Seq(sequence).reverse_complement())\n",
    "\n",
    "    nucleotide_pos= len(sequence)- (nucleotide_pos - 1)\n",
    "    nucleotide_pattern = re.compile(f\"(?P<before>.{{0,{nucleotide_pos-1}}})(?P<nucleotide>.)(?P<after>.*)\")\n",
    "    formatted_sequence = nucleotide_pattern.sub(\"\\g<before>\\033[1m\\g<nucleotide>\\033[0m\\g<after>\", revcomp_sequence)\n",
    "    print(f\"Here is reverse complement:                {formatted_sequence}\")\n",
    "if search_strand == \"sense\":\n",
    "    sequence=sequence\n",
    "else: \n",
    "    sequence=revcomp_sequence\n",
    "    \n",
    "# Find all PAM positions on the selected strand\n",
    "pam_positions = find_pam_positions(sequence if search_strand == \"sense\" else revcomp_sequence, search_strand)\n",
    "\n",
    "\n",
    "\n",
    "# Find all gRNAs \n",
    "sgRNA_list = []\n",
    "\n",
    "for pos in pam_positions:\n",
    "    if (search_strand == \"sense\" and 0 <= (pos - nucleotide_pos) < 20) or (search_strand == \"antisense\" and 0 <= (pos - nucleotide_pos) < 20):\n",
    "        if search_strand == \"sense\":\n",
    "            pPAM= pos + 3\n",
    "            gRNA_start = pos - 20\n",
    "            gRNA = sequence[gRNA_start:pos]\n",
    "            gPAM= sequence[gRNA_start:pPAM]\n",
    "        else:\n",
    "            pPAM= pos + 3\n",
    "            gRNA_start = pos - 20\n",
    "            gRNA = revcomp_sequence[gRNA_start:pos]\n",
    "            gPAM= sequence[gRNA_start:pPAM] \n",
    "            \n",
    "        cg_percent = calculate_cg_percent(gRNA)\n",
    "        \n",
    "        nucleotide_pos = nucleotide_pos  # Specify the nucleotide position of interest\n",
    "        end = pos + 1  #pos - 1\n",
    "        start = end - 20\n",
    "        pos_grna = 20 - (end - nucleotide_pos)\n",
    "        wt = gPAM  # Assign gRNA to wt\n",
    "\n",
    "        # Create a list of nucleotides in the gRNA with the specified replacements\n",
    "        off_list = [nt if i != pos_grna else 'G' if nt == 'A' else 'T' if nt == 'C' else nt for i, nt in enumerate(gPAM)]\n",
    "\n",
    "        # Join the modified nucleotides to form the off-gRNA sequence\n",
    "        off = ''.join(off_list)\n",
    "        # calculate cfd score \n",
    "        cfd_score= cfd_score_calculator(wt, off)\n",
    "\n",
    "\n",
    "        sorted_efficacies = find_sorted_efficacies(BE_profile, nucleotide_pos) \n",
    "        tm= primer3.calc_tm(gRNA)\n",
    "        Tm = round(tm, 2)\n",
    "        hp=primer3.calc_hairpin(gRNA)\n",
    "        Hairpin=hp.structure_found\n",
    "        \n",
    "        sgRNA_list.append([gRNA +' NGG', sorted_efficacies, cfd_score, cg_percent, Tm, Hairpin ])\n",
    "\n",
    "headers = [\"sgRNA Sequence\", \"Cell     Enzyme  Efficacy(%)\", \"CFD specificity score\", \"GC(%)\", 'Tm(°C)', \"Hairpin?\"]\n",
    "for row in sgRNA_list:\n",
    "    row[1] = tabulate(row[1], tablefmt=\"plain\")  # Format sorted efficacies as a table\n",
    "print(tabulate(sgRNA_list, headers=headers))\n",
    "\n",
    "print(\"----------------------------------------------------------------------------------------\")\n",
    "print(\"Here are restriction enzymes(if any) to test on wild and mutated sequences. \")\n",
    "\n",
    "#####Generate mutant sequence and find possible RFLP restriction enzyms\n",
    "generate_mutant_sequence(sequence, nucleotide_pos)\n",
    "\n",
    "\n",
    "print(\"----------------------------------------------------------------------------------------\")\n",
    "\n",
    "print('''\n",
    "\n",
    "Here are 10 primer pair combinations for ARMS PCR, suggested by Primer1(http://primer1.soton.ac.uk/primer1.html)\n",
    "''')\n",
    "run_doPrimer(org_sequence, org_nucleotide_pos, org_nucleotide, mut_nucleotide)\n",
    "\n",
    "\n",
    "print(\"----------------------------------------------------------------------------------------\")\n",
    "\n",
    "\n",
    "\n",
    "print ('''\n",
    "\n",
    "******************************************************Disclaimer:***********************************************\n",
    "This code is for educational purposes only. Use it at your own risk. The author of this code makes no warranties\n",
    "about the accuracy or suitability of this code for any purpose. By using this code, you agree to hold the author\n",
    "harmless from any damages that may arise from its use.\n",
    "''')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "30d66255",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: rs3 in /Users/alietemadi/opt/miniconda3/envs/BE3.7/lib/python3.7/site-packages (0.0.15)\n",
      "Requirement already satisfied: pip in /Users/alietemadi/opt/miniconda3/envs/BE3.7/lib/python3.7/site-packages (from rs3) (23.1.2)\n",
      "Requirement already satisfied: packaging in /Users/alietemadi/opt/miniconda3/envs/BE3.7/lib/python3.7/site-packages (from rs3) (23.1)\n",
      "Requirement already satisfied: joblib>=1.0.1 in /Users/alietemadi/opt/miniconda3/envs/BE3.7/lib/python3.7/site-packages (from rs3) (1.3.1)\n",
      "Requirement already satisfied: pandas>=1.0.0 in /Users/alietemadi/opt/miniconda3/envs/BE3.7/lib/python3.7/site-packages (from rs3) (1.3.5)\n",
      "Requirement already satisfied: lightgbm>=3.0.0 in /Users/alietemadi/opt/miniconda3/envs/BE3.7/lib/python3.7/site-packages (from rs3) (3.3.5)\n",
      "Requirement already satisfied: sglearn>=1.2.5 in /Users/alietemadi/opt/miniconda3/envs/BE3.7/lib/python3.7/site-packages (from rs3) (1.2.5)\n",
      "Requirement already satisfied: tqdm>=4.61.2 in /Users/alietemadi/opt/miniconda3/envs/BE3.7/lib/python3.7/site-packages (from rs3) (4.65.0)\n",
      "Requirement already satisfied: pyarrow>=4.0.1 in /Users/alietemadi/opt/miniconda3/envs/BE3.7/lib/python3.7/site-packages (from rs3) (12.0.1)\n",
      "Requirement already satisfied: biopython>=1.78 in /Users/alietemadi/opt/miniconda3/envs/BE3.7/lib/python3.7/site-packages (from rs3) (1.81)\n",
      "Requirement already satisfied: scikit-learn>=0.24.2 in /Users/alietemadi/opt/miniconda3/envs/BE3.7/lib/python3.7/site-packages (from rs3) (1.0.2)\n",
      "Requirement already satisfied: requests>=2.25.1 in /Users/alietemadi/opt/miniconda3/envs/BE3.7/lib/python3.7/site-packages (from rs3) (2.30.0)\n",
      "Requirement already satisfied: numpy in /Users/alietemadi/opt/miniconda3/envs/BE3.7/lib/python3.7/site-packages (from biopython>=1.78->rs3) (1.21.6)\n",
      "Requirement already satisfied: wheel in /Users/alietemadi/opt/miniconda3/envs/BE3.7/lib/python3.7/site-packages (from lightgbm>=3.0.0->rs3) (0.40.0)\n",
      "Requirement already satisfied: scipy in /Users/alietemadi/opt/miniconda3/envs/BE3.7/lib/python3.7/site-packages (from lightgbm>=3.0.0->rs3) (1.7.3)\n",
      "Requirement already satisfied: python-dateutil>=2.7.3 in /Users/alietemadi/opt/miniconda3/envs/BE3.7/lib/python3.7/site-packages (from pandas>=1.0.0->rs3) (2.8.2)\n",
      "Requirement already satisfied: pytz>=2017.3 in /Users/alietemadi/opt/miniconda3/envs/BE3.7/lib/python3.7/site-packages (from pandas>=1.0.0->rs3) (2023.3)\n",
      "Requirement already satisfied: charset-normalizer<4,>=2 in /Users/alietemadi/opt/miniconda3/envs/BE3.7/lib/python3.7/site-packages (from requests>=2.25.1->rs3) (3.1.0)\n",
      "Requirement already satisfied: idna<4,>=2.5 in /Users/alietemadi/opt/miniconda3/envs/BE3.7/lib/python3.7/site-packages (from requests>=2.25.1->rs3) (3.4)\n",
      "Requirement already satisfied: urllib3<3,>=1.21.1 in /Users/alietemadi/opt/miniconda3/envs/BE3.7/lib/python3.7/site-packages (from requests>=2.25.1->rs3) (1.26.7)\n",
      "Requirement already satisfied: certifi>=2017.4.17 in /Users/alietemadi/opt/miniconda3/envs/BE3.7/lib/python3.7/site-packages (from requests>=2.25.1->rs3) (2023.5.7)\n",
      "Requirement already satisfied: threadpoolctl>=2.0.0 in /Users/alietemadi/opt/miniconda3/envs/BE3.7/lib/python3.7/site-packages (from scikit-learn>=0.24.2->rs3) (3.1.0)\n",
      "Requirement already satisfied: seqfold>=0.7.7 in /Users/alietemadi/opt/miniconda3/envs/BE3.7/lib/python3.7/site-packages (from sglearn>=1.2.5->rs3) (0.7.17)\n",
      "Requirement already satisfied: six>=1.5 in /Users/alietemadi/opt/miniconda3/envs/BE3.7/lib/python3.7/site-packages (from python-dateutil>=2.7.3->pandas>=1.0.0->rs3) (1.16.0)\n"
     ]
    }
   ],
   "source": [
    "!pip install rs3\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "b503a7d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "from rs3.seq import predict_seq\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "cef88c21",
   "metadata": {},
   "outputs": [],
   "source": [
    "context_seqs = ['GACGAAAGCGACAACGCGTTCATCCGGGCA', 'AGAAAACACTAGCATCCCCACCCGCGGACT']\n",
    "predict_seq(context_seqs, sequence_tracr='Hsu2013')\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "52da93e2",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8554be76",
   "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}