SSDraw.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/ethanchen1301/SSDraw/blob/main/SSDraw.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "HEuP5CG7g9ak",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "%%capture\n",
        "#@title SSDraw\n",
        "#@markdown Installing the necessary dependencies can take 1-2 minutes.\n",
        "\n",
        "import os\n",
        "cwd=os.getcwd()\n",
        "print(cwd)\n",
        "os.chdir(\"/content/\")\n",
        "cwd=os.getcwd()\n",
        "print(cwd)\n",
        "\n",
        "# first install dssp\n",
        "#!sudo yes | apt-get install dssp\n",
        "!git clone https://github.com/ethanchen1301/SSDraw.git\n",
        "!tar -zxvf SSDraw/dssp-2.3.0.tar.gz\n",
        "%cd dssp-2.3.0\n",
        "!./autogen.sh\n",
        "!./configure\n",
        "!make\n",
        "!make mkdssp\n",
        "!sudo make install\n",
        "%cd ../\n",
        "!mkdssp --version\n",
        "\n",
        "# install necessary dependencies\n",
        "!pip install Bio\n",
        "!pip install matplotlib\n",
        "!pip install numpy\n",
        "\n",
        "import matplotlib.patches as mpatch\n",
        "import matplotlib.pyplot as plt\n",
        "from matplotlib.collections import PatchCollection\n",
        "import matplotlib.path as mpath\n",
        "import matplotlib\n",
        "import numpy as np\n",
        "import warnings\n",
        "from Bio import BiopythonDeprecationWarning\n",
        "warnings.simplefilter('ignore',BiopythonDeprecationWarning)\n",
        "from Bio import pairwise2\n",
        "import sys, os, string, math, re\n",
        "from matplotlib import rcParams\n",
        "import argparse\n",
        "from Bio.PDB import PDBParser\n",
        "from Bio.PDB.DSSP import DSSP\n",
        "from Bio.Align import substitution_matrices\n",
        "from Bio import AlignIO\n",
        "import matplotlib.colors as mcolors\n",
        "from matplotlib.colors import ListedColormap\n",
        "from Bio.SeqUtils import seq1\n",
        "import warnings, textwrap\n",
        "from Bio import BiopythonWarning\n",
        "from Bio.PDB.PDBExceptions import PDBConstructionWarning\n",
        "\n",
        "def read_r4s(inputfile):\n",
        "    # algorithm to convert raw scores to grades is taken from consurf:\n",
        "    # https://github.com/Rostlab/ConSurf\n",
        "    seq = \"\"\n",
        "    scores = []\n",
        "\n",
        "    pattern = r\"^\\s*?\\d+\\s+(\\w)\\s+(\\S+)\\s+\\[\\s*\\S+,\\s*\\S+\\]\\s+\\S+\\s+\\d+\\/\\d+\"\n",
        "    with open(inputfile) as READ:\n",
        "        for line in READ:\n",
        "            line = line.strip()\n",
        "            if re.match(pattern, line):\n",
        "                match = re.match(pattern, line)\n",
        "                seq+=match.group(1)\n",
        "                scores.append(eval(match.group(2)))\n",
        "\n",
        "    max_cons = min(scores)\n",
        "\n",
        "    # 9 steps from -|max_cons| to |max_cons|, midpoint is 0\n",
        "    ConsGradeUnity = max_cons / 4.5 * -1\n",
        "    if max_cons >= 0:\n",
        "        ConsGradeUnity = max_cons\n",
        "\n",
        "    grades = []\n",
        "    for score in scores:\n",
        "        grades.append(max(1, 9-int((score-max_cons)/ConsGradeUnity)))\n",
        "\n",
        "    return seq, grades\n",
        "\n",
        "def read_consurf_grad(inputfile):\n",
        "\n",
        "    grades = []\n",
        "    seq = \"\"\n",
        "    pattern = r\"^\\s*?\\d+\\s+(\\w)\\s+\\S+\\s+\\S+\\s+(\\d)\\S*\\s+-?\\d+.\\d+,\\s+-?\\d+.\\d+\\s+\\d,\\d\\s+\\d+\\/\\d+\\s+\\S+\"\n",
        "    with open(inputfile, \"r\") as f:\n",
        "        for line in f:\n",
        "            if re.match(pattern, line):\n",
        "                match = re.match(pattern,line)\n",
        "                seq += match.group(1)\n",
        "                grades.append(match.group(2))\n",
        "\n",
        "    return seq,grades\n",
        "\n",
        "def check_consurf_file(file):\n",
        "    consurf_pattern = r\"^\\s*?\\d+\\s+(\\w)\\s+\\S+\\s+\\S+\\s+(\\d)\\S*\\s+-?\\d+.\\d+,\\s+-?\\d+.\\d+\\s+\\d,\\d\\s+\\d+\\/\\d+\\s+\\S+\"\n",
        "    r4s_pattern = r\"^\\s*?\\d+\\s+(\\w)\\s+(\\S+)\\s+\\[\\s*\\S+,\\s*\\S+\\]\\s+\\S+\\s+\\d+\\/\\d+\"\n",
        "    with open(file, \"r\") as f:\n",
        "        for line in f:\n",
        "            if re.match(consurf_pattern, line):\n",
        "                return \"consurf\"\n",
        "            if re.match(r4s_pattern, line):\n",
        "                return \"r4s\"\n",
        "\n",
        "def gap_sequence(seq, extra_gaps):\n",
        "    # seq can be a list or a string, anything that can be indexed\n",
        "    # extra gaps is a list of length two [x,y],\n",
        "    # where x is the number of characters to remove from the beginning\n",
        "    # and y the number of characters to remove from the end\n",
        "    new_seq = seq\n",
        "    if extra_gaps[1] != 0:\n",
        "        new_seq = new_seq[:-extra_gaps[1]]\n",
        "    return new_seq[extra_gaps[0]:]\n",
        "\n",
        "def NormalizeData(data):\n",
        "    if np.min(data) == np.max(data):\n",
        "        warnings.warn(\"Warning: scores are the same for all residues\")\n",
        "        return data\n",
        "    return (data - np.min(data)) / (np.max(data) - np.min(data))\n",
        "\n",
        "def coords2path(coord_set1):\n",
        "\n",
        "    coords_f1 = []\n",
        "    instructions1 = []\n",
        "\n",
        "    for c in coord_set1:\n",
        "        for n in range(len(c)):\n",
        "            coords_f1.append(c[n])\n",
        "            if n == 0:\n",
        "                instructions1.append(1)\n",
        "            else:\n",
        "                instructions1.append(2)\n",
        "\n",
        "    return coords_f1, instructions1\n",
        "\n",
        "def build_loop(loop,idx,ssidx,loop_coords,linelen,nlines,prev_ss,next_ss,\n",
        "               z=1,clr='r',mat=0,size=75):\n",
        "\n",
        "    # if loop is smaller than 3 residues and has gaps on both sides, don't draw\n",
        "    if prev_ss == \"B\" and next_ss == \"B\" and loop[1]-loop[0] < 2:\n",
        "        return\n",
        "\n",
        "    i0 = loop[0]\n",
        "    if loop[0] != 0 and prev_ss != \"B\":\n",
        "        i0 = loop[0]-1\n",
        "    elif i0 == 0:\n",
        "        i0 = 0.06\n",
        "    else:\n",
        "        i0 = loop[0]\n",
        "    i1 = loop[1]+2\n",
        "    if loop[1] == linelen-1:\n",
        "        i1 += 2\n",
        "\n",
        "    o = 2\n",
        "    if idx == nlines-1:\n",
        "        o = 0\n",
        "    if next_ss == \"B\":\n",
        "        o = -1.5\n",
        "    if next_ss == None:\n",
        "        o=-4.1\n",
        "\n",
        "    rectangle = mpatch.Rectangle((i0/6.0,-0.25-5.5*idx+2.0*ssidx),\n",
        "                                 (i1-i0+o)/6.0,0.5,\n",
        "                                  fc=clr,ec='k',zorder=0)\n",
        "\n",
        "    plt.gca().add_patch(rectangle)\n",
        "\n",
        "    xy = rectangle.get_xy()\n",
        "    w = rectangle.get_width()\n",
        "    h = rectangle.get_height()\n",
        "    loop_coords.append(np.array([[xy[0],xy[1]],\n",
        "                                 [xy[0],xy[1]+h],\n",
        "                                 [xy[0]+w,xy[1]+h],\n",
        "                                 [xy[0]+w,xy[1]],\n",
        "                                 [xy[0],xy[1]]]))\n",
        "\n",
        "def build_strand(strand,idx,ssidx,strand_coords,next_ss,z=1,clr='r',\n",
        "                 imagemat=0,size=75):\n",
        "\n",
        "    delta = 0 if next_ss == None else 1\n",
        "\n",
        "    arrow=mpatch.FancyArrow(((strand[0]+delta-1)/6.0),-5.5*idx+2.0*ssidx,\n",
        "                            (strand[1]-strand[0]+1)/6.0,0,\n",
        "                            width=1.0,fc=clr,linewidth=0.5,ec='k',\n",
        "                            zorder=z,head_width=2.0,\n",
        "                            length_includes_head=True,head_length=2.0/6.0)\n",
        "\n",
        "    plt.gca().add_patch(arrow)\n",
        "\n",
        "    strand_coords.append(arrow.get_xy())\n",
        "\n",
        "\n",
        "def build_helix(helix,idx,ssidx,coord_set1, coord_set2, clr='r',size=37.5,\n",
        "                z=1,bkg=(0.195,0,0.051),imagemat=0):\n",
        "\n",
        "    i = helix\n",
        "    l = i[1]-i[0]+1\n",
        "    points = [[i[0]/6.0,-0.25-5.5*idx+2.0*ssidx],\n",
        "              [i[0]/6.0+1.0/6,0.75-5.5*idx+2.0*ssidx],\n",
        "              [i[0]/6.0+2.0/6,0.75-5.5*idx+2.0*ssidx],\n",
        "              [i[0]/6.0+1.0/6,-0.25-5.5*idx+2.0*ssidx]]\n",
        "    #hlx = plt.Polygon(points,fc=clr,ec='k',zorder=1,linewidth=2)\n",
        "    #coords= hlx.get_xy()\n",
        "    #coord_set2.append(coords)\n",
        "    coord_set2.append(points+[points[0]])\n",
        "\n",
        "    for j in range((l-2)-1):\n",
        "        if j % 2 == 0:\n",
        "            points = [[i[0]/6.0+(1.0+j)/6,0.75-5.5*idx+2.0*ssidx],\n",
        "                      [i[0]/6.0+(2.0+j)/6,0.75-5.5*idx+2.0*ssidx],\n",
        "                      [i[0]/6.0+(3.0+j)/6,-0.75-5.5*idx+2.0*ssidx],\n",
        "                      [i[0]/6.0+(2.0+j)/6,-0.75-5.5*idx+2.0*ssidx]]\n",
        "            coord_set1.append(points+[points[0]])\n",
        "            #hlx = mpatch.Polygon(points,fc=bkg,zorder=z)\n",
        "\n",
        "        else:\n",
        "            points = [[i[0]/6.0+(1.0+j)/6,-0.75-5.5*idx+2.0*ssidx],\n",
        "                      [i[0]/6.0+(2.0+j)/6,-0.75-5.5*idx+2.0*ssidx],\n",
        "                      [i[0]/6.0+(3.0+j)/6,0.75-5.5*idx+2.0*ssidx],\n",
        "                      [i[0]/6.0+(2.0+j)/6,0.75-5.5*idx+2.0*ssidx]]\n",
        "            coord_set2.append(points+[points[0]])\n",
        "            #hlx = mpatch.Polygon(points,fc=clr,zorder=0)\n",
        "\n",
        "\n",
        "    if (l-2-1)%2 == 1:\n",
        "\n",
        "        points = [[i[1]/6.0-1.0/6,-0.75-5.5*idx+2.0*ssidx],\\\n",
        "                  [i[1]/6.0,-0.75-5.5*idx+2.0*ssidx],\\\n",
        "                  [i[1]/6.0+1.0/6,0.25-5.5*idx+2.0*ssidx],\n",
        "                  [i[1]/6.0,0.25-5.5*idx+2.0*ssidx]]\n",
        "\n",
        "        coord_set2.append(points+[points[0]])\n",
        "        #hlx = mpatch.Polygon(points,fc=clr,zorder=0)\n",
        "\n",
        "    else:\n",
        "        points = [[i[1]/6.0-1.0/6,0.75-5.5*idx+2.0*ssidx],\n",
        "                  [i[1]/6.0,0.75-5.5*idx+2.0*ssidx],\\\n",
        "                  [i[1]/6.0+1.0/6,-0.25-5.5*idx+2.0*ssidx],\n",
        "                  [i[1]/6.0,-0.25-5.5*idx+2.0*ssidx]]\n",
        "        coord_set1.append(points+[points[0]])\n",
        "\n",
        "        #hlx = plt.Polygon(points,fc=bkg,zorder=10)\n",
        "\n",
        "def SS_breakdown(ss):\n",
        "\n",
        "    i = 0\n",
        "    curSS = ''\n",
        "    jstart = -1\n",
        "    jend = -1\n",
        "\n",
        "    strand = []\n",
        "    loop = []\n",
        "    helix = []\n",
        "    ssbreak = []\n",
        "\n",
        "    ss_order = []\n",
        "    ss_bounds = []\n",
        "\n",
        "    last_ss = ''\n",
        "\n",
        "    SS_equivalencies = {'H':['H'],\n",
        "                        '-':['-'],\n",
        "                        'S':[' ','S','C','T','G','I','P'],\n",
        "                        ' ':[' ','S','C','T','G','I','P'],\n",
        "                        'C':[' ','S','C','T','G','I','P'],\n",
        "                        'T':[' ','S','C','T','G','I','P'],\n",
        "                        'G':[' ','S','C','T','G','I','P'],\n",
        "                        'I':[' ','S','C','T','G','I','P'],\n",
        "                        'P':[' ','S','C','T','G','I','P'],\n",
        "                        'E':['E','B'],\n",
        "                        'B':['E','B']}\n",
        "\n",
        "    cur_SSDict = {'H':'helix',\n",
        "                  '-':'break',\n",
        "                  'E':'strand',\n",
        "                  'B':'strand'}\n",
        "\n",
        "    for i in range(len(ss)):\n",
        "\n",
        "        if i == 0:\n",
        "            curSS = SS_equivalencies[ss[i]]\n",
        "            jstart = i\n",
        "            if ss[i] in cur_SSDict.keys():\n",
        "                last_ss = cur_SSDict[ss[i]]\n",
        "            else:\n",
        "                last_ss = 'loop'\n",
        "            continue\n",
        "\n",
        "        if ss[i] in curSS:\n",
        "            jend = i\n",
        "\n",
        "        if ss[i] not in curSS or i == len(ss)-1:\n",
        "            if 'E' in curSS and jend-jstart+1 >= 3:\n",
        "                strand.append((jstart,jend))\n",
        "                ss_bounds.append((jstart,jend))\n",
        "                ss_order.append('E')\n",
        "                last_ss = 'strand'\n",
        "            elif 'H' in curSS and jend-jstart+1 >=4:\n",
        "                helix.append((jstart,jend))\n",
        "                ss_bounds.append((jstart,jend))\n",
        "                ss_order.append('H')\n",
        "                last_ss = 'helix'\n",
        "            elif ' ' in curSS and last_ss !='loop':\n",
        "                if jend < jstart:\n",
        "                    jend = jstart\n",
        "                loop.append((jstart,jend))\n",
        "                ss_bounds.append((jstart,jend))\n",
        "                ss_order.append('L')\n",
        "                last_ss = 'loop'\n",
        "            elif '-' in curSS:\n",
        "                if jend < jstart:\n",
        "                    jend = jstart\n",
        "                ssbreak.append((jstart,jend))\n",
        "                ss_bounds.append((jstart,jend))\n",
        "                ss_order.append('B')\n",
        "                last_ss = 'break'\n",
        "            elif last_ss == 'loop':\n",
        "                if jend < jstart:\n",
        "                    jend = jstart\n",
        "                if len(loop) > 0:\n",
        "                    jstart = loop[-1][0]\n",
        "                    loop = loop[0:-1]\n",
        "                    ss_bounds = ss_bounds[0:-1]\n",
        "                    ss_order = ss_order[0:-1]\n",
        "                loop.append((jstart,jend))\n",
        "                ss_bounds.append((jstart,jend))\n",
        "                ss_order.append('L')\n",
        "                last_ss = 'loop'\n",
        "            else:\n",
        "                if jend < jstart:\n",
        "                    jend = jstart\n",
        "                loop.append((jstart,jend))\n",
        "                ss_bounds.append((jstart,jend))\n",
        "                ss_order.append('L')\n",
        "                last_ss = 'loop'\n",
        "\n",
        "\n",
        "            jstart = i\n",
        "            curSS = SS_equivalencies[ss[i]]\n",
        "\n",
        "    return strand,loop,helix, ssbreak, ss_order, ss_bounds\n",
        "\n",
        "\n",
        "def updateSS(ss,seq,alignment,ref_align):\n",
        "\n",
        "    ss_u = ''\n",
        "\n",
        "    j = 0\n",
        "\n",
        "    for i in range(len(alignment)):\n",
        "        if alignment[i] == '-':\n",
        "            if ref_align[i] == '-':\n",
        "                ss_u += '-'\n",
        "            else:\n",
        "                ss_u += 'C'\n",
        "        else:\n",
        "            ss_u += ss[j]\n",
        "            j += 1\n",
        "\n",
        "    return ss_u\n",
        "\n",
        "\n",
        "def SS_align(alignment,ID,seq,ss,i_start,i_end):\n",
        "\n",
        "    a_seq = ''\n",
        "    seq_found = 0\n",
        "\n",
        "    for i in alignment:\n",
        "\n",
        "        if seq_found and i[0] == '>':\n",
        "            break\n",
        "\n",
        "        if i[0] == '>' and bool(re.search(ID.lower(), i.lower())):\n",
        "            seq_found = 1\n",
        "            continue\n",
        "\n",
        "        if seq_found and i[0] != '>':\n",
        "            a_seq += i\n",
        "\n",
        "    if i_end != 0:\n",
        "        i_end = len(a_seq)-i_end-1\n",
        "\n",
        "    a_seq = gap_sequence(a_seq, [i_start,i_end])\n",
        "\n",
        "\n",
        "    a = pairwise2.align.localxs(seq,a_seq,-1,-0.5)\n",
        "\n",
        "    # check if the dssp annotation has any extra residues not in the alignment\n",
        "    if a[0][1] != a_seq:\n",
        "        print(\"extra residues in pdb found\\n\")\n",
        "\n",
        "    # check how many gap marks are at the end and beginning of the alignment\n",
        "    # (a_seq) and compare to the amount found in a[0][1]\n",
        "    a_seq_gaps = [0,0]\n",
        "    new_aln_gaps = [0,0]\n",
        "    for i in range(len(a_seq)):\n",
        "        if a_seq[i] == \"-\":\n",
        "            a_seq_gaps[0] += 1\n",
        "        else:\n",
        "            break\n",
        "    for i in range(len(a_seq)-1, -1, -1):\n",
        "        if a_seq[i] == \"-\":\n",
        "            a_seq_gaps[1] += 1\n",
        "        else:\n",
        "            break\n",
        "\n",
        "    for i in range(len(a[0][1])):\n",
        "        if a[0][1][i] == \"-\":\n",
        "            new_aln_gaps[0] += 1\n",
        "        else:\n",
        "            break\n",
        "    for i in range(len(a[0][1])-1, -1, -1):\n",
        "        if a[0][1][i] == \"-\":\n",
        "            new_aln_gaps[1] += 1\n",
        "        else:\n",
        "            break\n",
        "\n",
        "    extra_gaps = [new_aln_gaps[0]-a_seq_gaps[0], new_aln_gaps[1]-a_seq_gaps[1]]\n",
        "\n",
        "    SS_updated = updateSS(ss,seq,a[0][0],a[0][1])\n",
        "\n",
        "    SS_updated_new = gap_sequence(SS_updated, extra_gaps)\n",
        "    a_new = gap_sequence(a[0][1], extra_gaps)\n",
        "\n",
        "    return SS_updated_new, a_new, extra_gaps, i_start, i_end\n",
        "\n",
        "\n",
        "def plot_coords(coords_all,mat,sz,CMAP):\n",
        "\n",
        "    for i,coords in enumerate(coords_all):\n",
        "\n",
        "        if not coords:\n",
        "            continue\n",
        "\n",
        "        coords_f1, instructions1 = coords2path(coords)\n",
        "\n",
        "        #If loop or bottom helix layer, zorder = 0\n",
        "        if i in [0,1]:\n",
        "            z = 0\n",
        "        else:\n",
        "            z = 10\n",
        "\n",
        "        path = mpath.Path(np.array(coords_f1),np.array(instructions1))\n",
        "        patch = mpatch.PathPatch(path, facecolor='none',zorder=z)\n",
        "        plt.gca().add_patch(patch)\n",
        "        im = plt.imshow(mat,extent=[0.0,sz,0.5,3.0],cmap=CMAP,interpolation='none',zorder=z)\n",
        "        im.set_clip_path(patch)\n",
        "\n",
        "\n",
        "def run_dssp(pdb_path, id, chain_id, dssp_exe='mkdssp'):\n",
        "\n",
        "    dssp_mode = \"dssp\"\n",
        "    ss_seq = \"\"\n",
        "    aa_seq = \"\"\n",
        "    #Suppress Biopython PDBParser warning\n",
        "    with warnings.catch_warnings(record=True) as w:\n",
        "        warnings.simplefilter(\"always\", PDBConstructionWarning)\n",
        "        p = PDBParser()\n",
        "        structure = p.get_structure(id, pdb_path)\n",
        "    model = structure[0]\n",
        "\n",
        "    dssp = DSSP(model, pdb_path, dssp=dssp_exe)\n",
        "\n",
        "\n",
        "    if dssp_mode == \"dssp\":\n",
        "        a_key = list(dssp.keys())\n",
        "        for key in a_key:\n",
        "            if key[0] == chain_id:\n",
        "                aa_seq+=dssp[key][1]\n",
        "                if dssp[key][2] == \"-\":\n",
        "                    ss_seq+=\"C\"\n",
        "                else:\n",
        "                    ss_seq+=dssp[key][2]\n",
        "\n",
        "    return [ss_seq,aa_seq]\n",
        "\n",
        "def convert2horiz(dssp_file, pdbseq):\n",
        "\n",
        "    ss_seq = \"\"\n",
        "    aa_seq = \"\"\n",
        "\n",
        "    with open(dssp_file, \"r\") as f:\n",
        "        lines = f.readlines()\n",
        "\n",
        "    dssp_ext = dssp_file.split(\".\")[-1]\n",
        "\n",
        "    if dssp_ext == \"horiz\":\n",
        "        ss_seq = lines[0].rstrip('\\n')\n",
        "        aa_seq = lines[1].rstrip('\\n')\n",
        "\n",
        "    elif dssp_ext == \"dssp\":\n",
        "        start_read = False\n",
        "        for line in lines:\n",
        "            if start_read:\n",
        "                if line[13] == \"!\":\n",
        "                    start_read = False\n",
        "                    continue\n",
        "                ss_seq += line[16]\n",
        "                aa_seq += line[13]\n",
        "            if line.split()[0] == \"#\":\n",
        "                start_read = True\n",
        "\n",
        "    else:\n",
        "        raise Exception(\"DSSP file extension not recognized: must be in .dssp or .horiz format\")\n",
        "\n",
        "    return [ss_seq,aa_seq]\n",
        "\n",
        "def score_column(msa_col, threshold=0):\n",
        "\n",
        "    blosum62 = substitution_matrices.load(\"BLOSUM62\")\n",
        "    # find consensus of the column\n",
        "    aa_count = {\"A\": 0, \"R\": 0, \"N\": 0, \"D\": 0, \"C\": 0, \"Q\": 0, \"E\": 0,\n",
        "                \"G\": 0, \"H\": 0, \"I\": 0, \"L\": 0, \"K\": 0, \"M\": 0, \"F\": 0,\n",
        "                \"P\": 0, \"S\": 0, \"T\": 0, \"W\": 0, \"Y\": 0, \"V\": 0}\n",
        "    for i in msa_col:\n",
        "        try:\n",
        "            aa_count[i] += 1\n",
        "        except:\n",
        "            pass\n",
        "    consensus_aa = max(zip(aa_count.values(), aa_count.keys()))[1]\n",
        "\n",
        "    conservation_count = 0\n",
        "    for i in msa_col:\n",
        "        if i in aa_count.keys():\n",
        "            if blosum62[consensus_aa][i] >= 0:\n",
        "                conservation_count += 1\n",
        "\n",
        "    return conservation_count/len(msa_col)\n",
        "\n",
        "def parse_color(args,seq_wgaps,pdbseq,bfactors,msa,extra_gaps):\n",
        "    CMAP = \"\"\n",
        "    if args.color in mcolors.BASE_COLORS.keys() or \\\n",
        "       args.color in mcolors.CSS4_COLORS.keys() or \\\n",
        "       args.color in mcolors.XKCD_COLORS.keys():\n",
        "        CMAP = ListedColormap([args.color])\n",
        "    elif args.color[0] == \"#\":\n",
        "        CMAP = ListedColormap([args.color])\n",
        "    if args.conservation_score or args.bfactor or args.scoring_file:\n",
        "        if len(args.color_map) == 1:\n",
        "            CMAP = args.color_map[0]\n",
        "        else:\n",
        "            CMAP = ListedColormap(args.color_map)\n",
        "\n",
        "    #bvals are to make the colormap; taken from input PDB\n",
        "    bvals = []\n",
        "\n",
        "    if args.mview:\n",
        "        mview_colors = {\"A\": 0, \"G\": 0, \"I\": 0, \"L\": 0, \"M\": 0, \"P\": 0, \"V\": 0,\n",
        "                        \"F\": 1, \"H\": 1, \"W\": 1, \"Y\": 1,\n",
        "                        \"K\": 2, \"R\": 2,\n",
        "                        \"D\": 3, \"E\": 3,\n",
        "                        \"S\": 4, \"T\": 4,\n",
        "                        \"N\": 5, \"Q\": 5,\n",
        "                        \"C\": 6}\n",
        "        mview_colors_hit = [0,0,0,0,0,0,0]\n",
        "\n",
        "        mview_color_map = [\"#33cc00\",\"#009900\",\"#cb0000\",\"#0133ff\",\"#0299fe\",\n",
        "                           \"#6601cc\",\"#ffff00\",\"#808080\"]\n",
        "\n",
        "        for i in range(len(seq_wgaps)):\n",
        "            try:\n",
        "                m = mview_colors[seq_wgaps[i]]\n",
        "                bvals.append(m)\n",
        "                mview_colors_hit[m]+=1\n",
        "            except:\n",
        "                bvals.append(7)\n",
        "\n",
        "        # remove colors of residues not in sequence\n",
        "        for i in range(len(mview_colors_hit)):\n",
        "            if mview_colors_hit[i] == 0:\n",
        "                mview_color_map.pop(i)\n",
        "                for j in range(len(bvals)):\n",
        "                    if bvals[j] > i:\n",
        "                        bvals[j] -= 1\n",
        "\n",
        "        CMAP = ListedColormap(mview_color_map)\n",
        "\n",
        "    elif args.consurf:\n",
        "        # read in a rate4site output file or consurf score file\n",
        "        # if rate4site output file, convert raw scores to 1-9 scores\n",
        "        consurf_color_map = [\"#10C8D1\", \"#8CFFFF\", \"#D7FFFF\", \"#EAFFFF\", \"#FFFFFF\",\n",
        "        \"#FCEDF4\", \"#FAC9DE\", \"#F07DAB\", \"#A02560\"]\n",
        "\n",
        "        scoring_seq = \"\"\n",
        "        bvals_tmp = []\n",
        "        consurf_mode = check_consurf_file(args.consurf)\n",
        "        if consurf_mode == \"consurf\":\n",
        "            scoring_seq,bvals_tmp = read_consurf_grad(args.consurf)\n",
        "        elif consurf_mode == \"r4s\":\n",
        "            scoring_seq,bvals_tmp = read_r4s(args.consurf)\n",
        "\n",
        "\n",
        "        # remove colors of residues not in sequence\n",
        "        for i in reversed(range(1,10)):\n",
        "            if i not in bvals_tmp:\n",
        "                consurf_color_map.pop(i-1)\n",
        "                for j in range(len(bvals_tmp)):\n",
        "                    if bvals_tmp[j] > i:\n",
        "                        bvals_tmp[j] -= 1\n",
        "\n",
        "        CMAP = ListedColormap(consurf_color_map)\n",
        "\n",
        "        score_align = pairwise2.align.localxs(pdbseq,scoring_seq,-1,-0.5)\n",
        "\n",
        "        j = 0\n",
        "        for i in range(len(score_align[0][1])):\n",
        "            if score_align[0][0][i] != \"-\":\n",
        "                if score_align[0][1][i] != \"-\":\n",
        "                    bvals.append(bvals_tmp[j])\n",
        "                    j+=1\n",
        "                else:\n",
        "                    bvals.append(min(bvals_tmp))\n",
        "\n",
        "    elif args.scoring_file: # use custom scoring by residue\n",
        "        # read in scoring file\n",
        "        bvals_tmp = []\n",
        "        scoring_seq = \"\"\n",
        "        with open(args.scoring_file, \"r\") as g:\n",
        "            lines = g.readlines()\n",
        "\n",
        "        if lines[0] == '\\t Amino Acid Conservation Scores\\n': # read in consurf file\n",
        "            for line in lines:\n",
        "                if len(line.split()) > 0:\n",
        "                    if line.split()[0].isdigit():\n",
        "                        scoring_seq += line.split()[1]\n",
        "                        bvals_tmp.append(float(line.split()[4]))\n",
        "\n",
        "        else: # read in custom scoring file\n",
        "            for line in lines:\n",
        "                scoring_seq += line.split()[0]\n",
        "                bvals_tmp.append(float(line.split()[1]))\n",
        "\n",
        "        score_align = pairwise2.align.localxs(pdbseq,scoring_seq,-1,-0.5)\n",
        "\n",
        "        j = 0\n",
        "        for i in range(len(score_align[0][1])):\n",
        "            if score_align[0][0][i] != \"-\":\n",
        "                if score_align[0][1][i] != \"-\":\n",
        "                    bvals.append(bvals_tmp[j])\n",
        "                    j+=1\n",
        "                else:\n",
        "                    bvals.append(min(bvals_tmp))\n",
        "\n",
        "\n",
        "    elif args.bfactor:  # score by bfactor\n",
        "        bvals = [b for b in bfactors]\n",
        "\n",
        "    elif args.conservation_score: # score by conservation score\n",
        "        bvals = []\n",
        "        for i in range(len(msa[0])):\n",
        "            bvals.append(score_column([msa[j][i] for j in range(len(msa))]))\n",
        "\n",
        "    else: # solid color\n",
        "        bvals = [i for i in range(len(msa[0]))]\n",
        "\n",
        "    if len(bvals) == len(pdbseq):\n",
        "        # remove extra residues\n",
        "        pdbseq = gap_sequence(pdbseq, extra_gaps)\n",
        "        bvals = gap_sequence(bvals, extra_gaps)\n",
        "\n",
        "        j = 0\n",
        "        bvalsf = []\n",
        "        bvals_align = pairwise2.align.localxs(seq_wgaps,pdbseq,-1,-0.5)\n",
        "        if bvals_align[0][0] != seq_wgaps:\n",
        "            print(\"Error in alignment or pdb sequence\")\n",
        "            sys.exit()\n",
        "        for i in range(len(seq_wgaps)):\n",
        "            if seq_wgaps[i] == \"-\" or bvals_align[0][1][i] == \"-\":\n",
        "                if j >= len(bvals):\n",
        "                    bvalsf.append(bvals[j-1])\n",
        "                else:\n",
        "                    bvalsf.append(bvals[j])\n",
        "            elif seq_wgaps[i] == bvals_align[0][1][i]:\n",
        "                bvalsf.append(bvals[j])\n",
        "                j+=1\n",
        "\n",
        "        bvals = bvalsf\n",
        "\n",
        "    return CMAP, bvals\n",
        "\n",
        "def read_pdb(id,args):\n",
        "    pdbseq = ''\n",
        "    #Suppress Biopython PDBParser warning\n",
        "    with warnings.catch_warnings(record=True) as w:\n",
        "        warnings.simplefilter(\"always\", PDBConstructionWarning)\n",
        "        p = PDBParser()\n",
        "        bfactors = []\n",
        "        structure = p.get_structure(id, args.pdb)\n",
        "        model = structure[0]\n",
        "        for chain in model:\n",
        "            if chain.get_id() == args.chain_id:\n",
        "                for residue in chain:\n",
        "                    for atom in residue:\n",
        "                        if atom.name == \"CA\":\n",
        "                            bfactors.append(atom.bfactor)\n",
        "                            pdbseq += seq1(residue.get_resname())\n",
        "                            break\n",
        "    sys.stdout = sys.__stdout__\n",
        "\n",
        "    return bfactors, pdbseq\n",
        "\n",
        "description='''\n",
        "----------------\n",
        "\n",
        "SSDraw is a program that generates publication-quality protein secondary structure diagrams from three-dimensional protein structures. To depict relationships between secondary structure and other protein features, diagrams can be colored by conservation score, B-factor, or custom scoring.\n",
        "\n",
        "SSDraw also has a colab notebook available at https://github.com/ethanchen1301/SSDraw/blob/main/SSDraw.ipynb\n",
        "\n",
        "\n",
        "\n",
        "Installation:\n",
        "\n",
        "SSDraw requires the biopython module:\n",
        "        pip install biopython\n",
        "\n",
        "SSDraw also requires the DSSP program to be installed in order to generate secondary structure annotations.\n",
        "        sudo apt-get install dssp\n",
        "\n",
        "Alternatively, you can install DSSP either through conda (conda install -c salilab dssp), or you can follow the instructions on their github page to make a local installation:\n",
        "https://github.com/cmbi/dssp.\n",
        "\n",
        "\n",
        "\n",
        "Instructions:\n",
        "SSDraw requires 4 arguments:\n",
        "1. --fasta: the file name sequence or alignment file in fasta format.\n",
        "2. --name: the id of the sequence in the fasta file corresponding to your protein of interest.\n",
        "3. --pdb: the file name of the pdb file of your protein\n",
        "4. --output: the output file name to use\n",
        "\n",
        "Example 1:\n",
        "    python3 ../SSDraw.py --fasta 1ndd.fasta --name 1ndd --pdb 1ndd.pdb --output 1ndd_out\n",
        "\n",
        "Coloring options:\n",
        "SSDraw uses a gradient to color each position in the alignment by a certain score. The user can choose which scoring system to use, and they can also choose which colormap.\n",
        "\n",
        "Scoring:\n",
        "-conservation_score: score each position in the alignment by conservation score.\n",
        "-bfactor: score each residue in the pdb by B-factor\n",
        "-scoring_file: score each residue by a custom scoring file prepared by the user\n",
        "-mview: color each residue by the mview coloring system\n",
        "\n",
        "Example 2: Score by conservation\n",
        "    python3 ../SSDraw.py --fasta aligned.fasta --name 1ndd --pdb 1ndd.pdb --output 1ndd_conservation -conservation_score\n",
        "\n",
        "Example 3: Score by B-factor\n",
        "    python3 ../SSDraw.py --fasta 1ndd.fasta --name 1ndd --pdb 1ndd.pdb --output 1ndd_bfactor -bfactor\n",
        "\n",
        "Choosing a colormap:\n",
        "The default colormap for SSDraw is inferno. The user can select one of the matplotlib library color maps or simply list a set of colors they'd like to use with the --color_map option. Alternatively, the user can select a single color with the --color option and SSDraw will use that color on the whole image.\n",
        "\n",
        "Example 4: Custom scoring file with custom color map\n",
        "    python3 ../SSDraw.py --fasta 2kdl.fasta --name 2kdl --pdb 2kdl.pdb --output 2kdl_out --scoring_file 2kdl_scoring.txt --color_map black cyan\n",
        "\n",
        "DSSP files:\n",
        "Normally, SSDraw will generate a DSSP annotation from the PDB file, but if you have a DSSP file you would like to use, you can upload it and input the file name in Options.\n",
        "\n",
        "--start and --end:\n",
        "If you want SSDraw to draw only a portion of your alignment, you can specify the start and/or end points using the --start and --end options respectively. The argument for these options correspond to the index of the alignment position, not to the residue position numbers.\n",
        "\n",
        "Example 5: Choose subregion of alignment to run SSDraw on\n",
        "    python3 ../SSDraw.py --fasta aligned.fasta --name 1ndd --pdb 1ndd.pdb --output 1ndd_conservation_cropped -conservation_score --start 80 --end 132\n",
        "\n",
        "Running on multiple pdbs:\n",
        "In order to rapidly generate multiple images with SSDraw, we recommend writing shell scripts comprised of commands like those shown in the above examples. For examples of such shell scripts, see one of the shell scripts in /figures/.\n",
        "----------------\n",
        "\n",
        "'''"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "&#8592;\n",
        "\n",
        "&#8592;\n",
        "\n",
        "You will need to upload the sequence and pdb files to this session.\n",
        "Click on the Files icon to the left then choose 'upload to session storage'\n",
        "\n",
        "Once the files are uploaded type the file names below.\n",
        "\n",
        "&#8592;\n",
        "\n",
        "&#8592;"
      ],
      "metadata": {
        "id": "2kRtLMtIo7HF"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "EoaVx2yihGke",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "#%%capture\n",
        "\n",
        "## allow user to upload from their computer\n",
        "\n",
        "from google.colab import files\n",
        "#@title Define variables\n",
        "#@markdown ##Required:\n",
        "#@markdown Please input the following arguments:\n",
        "fasta_file = \"\" #@param {type:\"string\"}\n",
        "#@markdown  - The filename of the input sequence/alignment file (must be in fasta format).\n",
        "name = ''  #@param {type:\"string\"}\n",
        "#@markdown - The fasta identifier of target protein (what follows the \">\" character in the fasta file).\n",
        "pdb_file = \"\" #@param {type:\"string\"}\n",
        "#@markdown  - The filename of the input pdb file.\n",
        "output = \"\" #@param {type:\"string\"}\n",
        "#@markdown - The desired name of the output file (default image type is png).\n",
        "chain_id = 'A'  #@param {type:\"string\"}\n",
        "#@markdown - The chain id in the pdb file (default \"A\")\n",
        "\n",
        "#@markdown For more information please see [instructions](#Instructions) at the bottom of the page.\n",
        "\n",
        "salign = open(fasta_file).read().splitlines() # every line from the alignment\n",
        "\n",
        "#@markdown ##Options:\n",
        "#@markdown These are optional parameters you can use to further customize the generated image. You should not need to modify any of these to get the code to run.\n",
        "coloring = \"solid\" #@param [\"solid\", \"gradient\", \"mview\"]\n",
        "#@markdown - Choose whether to color the image by a gradient or with a solid coloring.\n",
        "scoring = \"by conservation score\" #@param [\"by conservation score\",\"by bfactor\", \"consurf\",\"custom\"]\n",
        "#@markdown - If coloring by gradient, the image will be colored according to the chosen scoring system.\n",
        "consurf = \"\" #@param {type: \"string\"}\n",
        "#@markdown - If you want to color by consurf grades, type in the name of your consurf file here\n",
        "scoring_file = \"\" #@param {type: \"string\"}\n",
        "#@markdown - If custom scoring was selected, please input the filename of the scoring file.\n",
        "gradient_color_map = \"inferno\" #@param {type: \"string\"}\n",
        "#@markdown - Choose which gradient color map to use (default is inferno). You can use any built-in colormap from matplotlib. You can also input a list of space-delimited colors (eg. \"black cyan red\").\n",
        "solid_color = \"white\" #@param {type:\"string\"}\n",
        "#@markdown - If coloring by solid, this is the color that will be used.\n",
        "dssp_file =\"\" #@param {type:\"string\"}\n",
        "#@markdown - If desired, one can choose to input their own secondary structure annotation for SSDraw to use. File must be in .dssp or .horiz format. Using this option is not recommended for most cases, as SSDraw already computes the secondary structure with DSSP from the given PDB file.\n",
        "ticks=0 #@param [0,5,10,20]\n",
        "#@markdown - Add ticks at every nth position to the output image. If this option is set to zero, ticks will not be added to the image.\n",
        "output_file_type = \"png\" #@param [\"png\", \"ps\", \"eps\", \"tif\", \"svg\"] {type:\"string\"}\n",
        "image_dpi = 600 #@param {type:\"raw\"}\n",
        "\n",
        "\n",
        "id = name\n",
        "color = solid_color\n",
        "color_map = gradient_color_map\n",
        "\n",
        "f = run_dssp(pdb_file, id, chain_id)\n",
        "\n",
        "# custom arg parse for colab notebook\n",
        "class ArgParser():\n",
        "    def __init__(self, fasta, name, pdb, output, dssp=None,\n",
        "                 chain_id = \"A\", color_map=[\"inferno\"], scoring_file=None,\n",
        "                 color=\"white\", conservation_score=False, output_file_type=\"png\",\n",
        "                 bfactor=False,mview=False,dpi=600,start=0,end=0,ticks=0):\n",
        "        self.fasta = fasta\n",
        "        self.name = name\n",
        "        self.pdb = pdb\n",
        "        self.output = output\n",
        "        self.dssp = dssp\n",
        "        self.chain_id = chain_id\n",
        "        self.color_map = color_map\n",
        "        self.scoring_file = scoring_file\n",
        "        self.color = color\n",
        "        self.conservation_score = conservation_score\n",
        "        self.output_file_type = output_file_type\n",
        "        self.bfactor = bfactor\n",
        "        self.mview = mview\n",
        "        self.dpi = dpi\n",
        "        self.start = start\n",
        "        self.end = end\n",
        "        self.ticks = int(ticks)\n",
        "        self.consurf = \"\"\n",
        "\n",
        "# parse arguments\n",
        "args = ArgParser(fasta_file,name,pdb_file,output,dssp=dssp_file,\n",
        "                   chain_id=chain_id,color_map=gradient_color_map,\n",
        "                   scoring_file=None, color=solid_color,\n",
        "                 conservation_score=False, output_file_type=output_file_type,\n",
        "                 bfactor=False,mview=False,dpi=600,start=0,end=0,ticks=ticks)\n",
        "\n",
        "# parse coloring and scoring\n",
        "if coloring != \"solid\":\n",
        "    args.color_map = gradient_color_map.split()\n",
        "    if coloring == \"mview\":\n",
        "        args.mview = True\n",
        "    elif scoring == \"consurf\":\n",
        "        args.consurf = consurf\n",
        "    elif scoring == \"by conservation score\":\n",
        "        args.conservation_score = True\n",
        "    elif scoring == \"by bfactor\":\n",
        "        args.bfactor = True\n",
        "    elif scoring == \"custom\":\n",
        "        args.scoring_file = scoring_file"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "id3iqycepcxe",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "#@title Choose subregion of alignment to run SSDraw on\n",
        "#@markdown If you want to run SSDraw on the entire alignment, run this cell, but leave these sliders unchanged.\n",
        "\n",
        "#@markdown IMPORTANT: These are the numberings based on the alignment, not on the numbering in the pdb file.\n",
        "## Change numbering to the numbering in the pdb\n",
        "\n",
        "# alignment subsection\n",
        "from ipywidgets import interact, interactive, fixed, interact_manual\n",
        "import ipywidgets as widgets\n",
        "from IPython.display import display\n",
        "\n",
        "display_seq = \"\"\n",
        "display_seq_found = False\n",
        "for i in salign:\n",
        "        if display_seq_found and i[0] == '>':\n",
        "            break\n",
        "        if i[0] == '>' and bool(re.search(id.lower(), i.lower())):\n",
        "            display_seq_found = True\n",
        "            continue\n",
        "\n",
        "        if display_seq_found and i[0] != '>':\n",
        "            display_seq += i\n",
        "\n",
        "def aln_sub(start=0,end=len(display_seq)):\n",
        "  display(display_seq[start:end+1])\n",
        "  return [start,end]\n",
        "\n",
        "w = interactive(aln_sub, start=(0,len(display_seq)), end=(0,len(display_seq)))\n",
        "display(w)"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "#@title Run SSDraw\n",
        "\n",
        "# parse start and end\n",
        "start,end = w.result\n",
        "args.start = start\n",
        "if end == len(display_seq):\n",
        "    args.end = 0\n",
        "else:\n",
        "    args.end = end\n",
        "\n",
        "def get_args():\n",
        "\n",
        "    parser = argparse.ArgumentParser(\n",
        "        formatter_class=argparse.RawDescriptionHelpFormatter,\n",
        "        description = description,\n",
        "        epilog=\"\")\n",
        "    parser.add_argument(\"-f\", \"--fasta\", help=\"(required) sequence/alignment file in fasta format\")\n",
        "    parser.add_argument(\"-p\", \"--pdb\", help=\"(required) pdb file\")\n",
        "    parser.add_argument(\"-n\", \"--name\", help=\"(required) id of the protein in the alignment file\")\n",
        "    parser.add_argument(\"-o\", \"--output\", help=\"(required) name for output file\")\n",
        "    parser.add_argument(\"--dssp\", default=None, help=\"secondary structure annotation in DSSP or .horiz format. If this option is not provided, SSDraw will compute secondary structure from the given PDB file with DSSP.\")\n",
        "    parser.add_argument(\"--chain_id\", default=\"A\", help=\"chain id to use in pdb. Defaults to chain A.\")\n",
        "    parser.add_argument(\"--color_map\", default=[\"inferno\"], nargs=\"*\", help=\"color map to use for heat map\")\n",
        "    parser.add_argument(\"--scoring_file\",default=None,help=\"custom scoring file for alignment\")\n",
        "    parser.add_argument(\"--color\", default=\"white\", help=\"color for the image. Can be a color name (eg. white, black, green), or a hex code\")\n",
        "    parser.add_argument(\"-conservation_score\", action='store_true', help=\"score alignment by conservation score\")\n",
        "    parser.add_argument(\"--output_file_type\", default=\"png\")\n",
        "    parser.add_argument(\"-bfactor\", action='store_true', help=\"score by B-factor\")\n",
        "    parser.add_argument(\"-mview\", action=\"store_true\", help=\"color by mview color map\")\n",
        "    parser.add_argument(\"--dpi\", default=600, type=int, help=\"dpi to use for final plot\")\n",
        "    parser.add_argument(\"--start\", default=0, type=int)\n",
        "    parser.add_argument(\"--end\", default=0,type=int)\n",
        "    parser.add_argument(\"--consurf\", default=\"\", help=\"consurf or rate4site file to color image with. If rate4site file is given, SSDraw will convert raw scores to grades.\")\n",
        "\n",
        "    args = parser.parse_args()\n",
        "\n",
        "    return args, parser\n",
        "\n",
        "def initialize(args):\n",
        "    #args, parser = get_args()\n",
        "\n",
        "    '''if not args.fasta or not args.pdb or not args.output or not args.name:\n",
        "        parser.print_help(sys.stderr)\n",
        "        sys.exit(1)'''\n",
        "\n",
        "    # preface run\n",
        "    id = args.name\n",
        "    chain_id = args.chain_id\n",
        "    print(\"\\nRunning for: \"+id)\n",
        "\n",
        "    # read in amino acid sequence from PDB\n",
        "    bfactors, pdbseq = read_pdb(id,args)\n",
        "\n",
        "    if args.dssp:\n",
        "        # get secondary structure from pre-existing DSSP annotation\n",
        "        f = convert2horiz(args.dssp,pdbseq)\n",
        "    else:\n",
        "        # run the dssp executable\n",
        "        f = run_dssp(args.pdb, id, chain_id)\n",
        "\n",
        "    nlines = 1\n",
        "    salign = open(args.fasta).read().splitlines()\n",
        "\n",
        "    if args.start > args.end:\n",
        "        raise Exception(\"--start cannot be greater than --end\")\n",
        "\n",
        "    #####Align secondary structure to match input sequence alignment\n",
        "    ss_wgaps,seq_wgaps,extra_gaps,i_start,i_end = SS_align(salign,args.name,f[1],f[0],args.start,args.end)\n",
        "\n",
        "    #Break down secondary structure classifications in to continuous\n",
        "    #chunks of helix, strand, and coil\n",
        "    strand,loop,helix,ss_break,ss_order,ss_bounds = SS_breakdown(ss_wgaps)\n",
        "\n",
        "    msa = [gap_sequence(a, [i_start, i_end]) for a in AlignIO.read(open(args.fasta), \"fasta\")]\n",
        "\n",
        "    return args,pdbseq,bfactors,msa,ss_wgaps,seq_wgaps,extra_gaps,i_start,i_end,strand,loop,helix,ss_break,ss_order,ss_bounds\n",
        "\n",
        "def SSDraw(args):\n",
        "\n",
        "    args,pdbseq,bfactors,msa,ss_wgaps,seq_wgaps,extra_gaps,i_start,i_end,strand,loop,helix,ss_break,ss_order,ss_bounds=initialize(args)\n",
        "    nlines = 1\n",
        "\n",
        "    #Parse color and scoring args\n",
        "    CMAP, bvals = parse_color(args,seq_wgaps,pdbseq,bfactors,msa,extra_gaps)\n",
        "\n",
        "    mat = np.tile(NormalizeData(bvals), (100,1))\n",
        "\n",
        "    #set figure parameters\n",
        "    sz = 0\n",
        "    c = 'none'\n",
        "    bc = 'none'\n",
        "\n",
        "    #set sizes of SS chunks\n",
        "    ss_prev = 0\n",
        "    for i in range(len(ss_order)):\n",
        "\n",
        "        if ss_order[i] == 'H':\n",
        "            ss_prev = ss_bounds[i][1]/6.0+1/6.0\n",
        "        else:\n",
        "            ss_prev = ss_bounds[i][1]/6.0\n",
        "\n",
        "    if ss_order[-1] == 'H':\n",
        "        sz = ss_bounds[-1][1]/6.0+1/6.0\n",
        "    elif ss_order[-1] in ['E','B']:\n",
        "        sz = ss_bounds[-1][1]/6.0\n",
        "    elif ss_order[-1] == 'L':\n",
        "        sz = (ss_bounds[-1][1])/6.0\n",
        "\n",
        "    #Plot secondary structure chunks\n",
        "    strand_coords = []\n",
        "    loop_coords = []\n",
        "    helix_coords1 = []\n",
        "    helix_coords2 = []\n",
        "\n",
        "    fig, ax = plt.subplots(ncols=1, figsize=(25,2+1.5*(nlines-1)))\n",
        "\n",
        "    for i in range(len(ss_order)):\n",
        "        prev_ss = None\n",
        "        next_ss = None\n",
        "        if i != 0:\n",
        "            prev_ss = ss_order[i-1]\n",
        "        if i != len(ss_order)-1:\n",
        "            next_ss = ss_order[i+1]\n",
        "\n",
        "        if ss_order[i] == 'L':\n",
        "            build_loop(ss_bounds[i],0,1,loop_coords,len(ss_wgaps),1,prev_ss,next_ss,z=0,clr=c,mat=mat,size=sz)\n",
        "        elif ss_order[i] == 'H':\n",
        "            build_helix(ss_bounds[i],0,1,helix_coords1,helix_coords2,z=i,clr=c,bkg=bc,imagemat=mat,size=sz)\n",
        "        elif ss_order[i] == 'E':\n",
        "            build_strand(ss_bounds[i],0,1,strand_coords,next_ss,z=i,clr=c,imagemat=mat,size=sz)\n",
        "\n",
        "\n",
        "    plot_coords([loop_coords,helix_coords2,strand_coords,helix_coords1],mat,sz,CMAP)\n",
        "\n",
        "    plt.ylim([0.5,3])\n",
        "\n",
        "    # remove spines and yticks\n",
        "    ax.spines['top'].set_visible(False)\n",
        "    ax.spines['right'].set_visible(False)\n",
        "    ax.spines['bottom'].set_visible(False)\n",
        "    ax.spines['left'].set_visible(False)\n",
        "    ax.get_yaxis().set_ticks([])\n",
        "\n",
        "    if args.ticks == 0:\n",
        "        ax.get_xaxis().set_ticks([])\n",
        "    else:\n",
        "        res_x = 0.1646\n",
        "        ticks = []\n",
        "        labels = []\n",
        "        i = 0\n",
        "        label_i = 1\n",
        "        while label_i <= len(bvals):\n",
        "            ticks.append(i)\n",
        "            labels.append(str(label_i))\n",
        "            i+=res_x*args.ticks\n",
        "            label_i+=args.ticks\n",
        "        ax.get_xaxis().set_ticks(ticks,labels=labels)\n",
        "        ax.xaxis.set_ticks_position('top')\n",
        "\n",
        "    ax.set_aspect(0.5)\n",
        "\n",
        "    print(\"Saving output to output/{:}.{:}...\".format(args.output, args.output_file_type))\n",
        "    if not os.path.exists(\"output\"):\n",
        "        os.makedirs(\"output\")\n",
        "    plt.savefig(\"output/\"+args.output+\".\"+args.output_file_type,bbox_inches='tight',dpi=args.dpi,transparent=True)\n",
        "\n",
        "\n",
        "SSDraw(args)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 106
        },
        "cellView": "form",
        "id": "apdnrtPLXcId",
        "outputId": "6fadee26-5a9a-45ea-9918-80364c0dfd89"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 2500x200 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAB6UAAACrCAYAAAAuGkkpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABL2UlEQVR4nO3deVzV5533/zc7iigiYFQURQQXQEA4C6DTtDOTbknaNCZp2klmptPOb7qkd5rMnaTtTNvMNN23tE2m7dztdDr3dLGZNLmbpUnTNkr8soiAC4oCIhB3UUGU9ZzfHwTUsJztOh4OvJ6PR/7ge67zOR8BJZzP93pfEW632y0AAAAAAAAAAAAAAIIgMtQNAAAAAAAAAAAAAABmLobSAAAAAAAAAAAAAICgYSgNAAAAAAAAAAAAAAgahtIAAAAAAAAAAAAAgKBhKA0AAAAAAAAAAAAACBqG0gAAAAAAAAAAAACAoGEoDQAAAAAAAAAAAAAIGobSAAAAAAAAAAAAAICgYSgNAAAAAAAAAAAAAAgahtIAAMxSLS0t+rvim1S4IltdXV2hbgeAn27MKtCt652qr68PdSsAAGCW6O/v163rnXrL6nwdOXIk1O0AAIBp6vnnn9fnPnK/CnLzNTg4GOp2EGIMpQEAmMXiomP0jb/4W33w7g+EuhUAAfj62z+ke7a+n1/wAADAdRMh6Tvv/Hvd8Z73yeVyhbodAAAwTaUlp+jBu/5OjzzySKhbQYgxlAYAYJbbtCxLOT3x2rZtW6hbAeCntHlJuj//7XrooYdC3QoAAJhFVi5crDuWFulrX/taqFsBAADT2J1/ebNam5rV0NAQ6lYQQgylAQCAPn3jnXrsgc8Q4w2Esds2lKnlRYsYbwAAcF19pOTteuGH/5cYbwAAMKUnvvEd3XvPPaS8zWIMpQEAgGIjo/XNv/wQMd5AmPvmuz5MjDcAALjuvvvufyDGGwAATCk1JUX/+In7ifGexRhKAwAASVLR0tXEeANhLmXOfGK8AQDAdbdiQSox3gAAwKO7bttKjPcsxlAaAACMIcYbCH/EeAMAgFAgxhsAAHiDGO/Zi6E0AAAYQ4w3MDMQ4w0AAEKBGG8AAOAJMd6zF0NpAABwDWK8gfBHjDcAAAgFYrwBAIA3iPGenRhKAwCAcYjxBsIfMd4AACAUiPEGAADeIMZ79mEoDQAAxiHGG5gZiPEGAACh8N13/4PuuH0rMd4AAGBSxHjPPgylAQDAhIjxBsIfMd4AACAUVixI1QffejMx3gAAYErEeM8uDKUBAMCkiPEGwh8x3gAAIBQ+uvUe/e6pZ4nxBgAAUyLGe/ZgKA0AACZFjDcwMxDjDQAAQuEHn/kyMd4AAGBKxHjPHtGhbgAAAExvRUtXK+dApbZt26atW7eGuh0Afrg6xvub3/xmqNtBiN17zz068IcqxUZHebU+MjnR45qTF86qqaU50NYkSS+88II+/N67NTciVpERER7Xp63N8LjmZNcZPfnjH+mtb31rwP253W7l5uYqKSlJ0dFmfqV2uVxKS0vTU089ZaTeAw88oD+98gfNS5hnpJ4ktXW0q/VIq6KivPu+mUpdXZ3e/76tSk5MUrSHeq6uHq9qdl3u1UNffVT33ntvwP0BMGflkvSxGG+OEwEAAJO567ateure36ihoUEbN24MdTsIEobSAADAo0/feKf+8oHP6G1ve5uSk5ND3Q4AP9y2oUy/+fU3VH9PvQoKCkLdDkJkYGBAf3ruJf389n/UhsWeh7mSNPfdmzyu+dvPP6DTp08rNTU10Bb13HPPKW/OMt2VWuLV+rv/7Z89rnml+jX96U9/MjKUbmtr08WLF7V///6Aa13N5N/L5557Tt/+4lf1Fze+zVjNhz//WR04cEC5ubkB1/rjH/+oTTm5+o8vfMvj2ku/rfWq5oFTHXrilVcYSgPT0Ee33qObPv5B3XHHHVq1alWo2wEAANPUE9/4jt6+9RbV7NqlmJiYULeDICC+GwAAeESMNzAzEOONl19+WVuWr1d15yGjdW25haqsrDRSq/VwszoHzhmpNcq2YaOampqM1LIsS2vWrFFHR4eReqNWr16t3t7egOsMDAwoeUGSKmuqDXR1hb3YJsuyjNTau3evjp0+ZaTWqJzUZTrR/rrRmgDMIcYbAAB4Qoz3zMdQGgAAeKVo6Wrl9MRr27ZtoW4FgJ+ujvHG7PTrbdv0sdKbtet1s0Pp8sISIwNLl8ulY7sPKVqRGnabG1wkJszTmY7jRmpZlqV3vetd2rFjh5F6o+x2u2pqagKu09DQoLeUbdHuhjoDXV1RZncYG0ofPdyi9MVLdK77gpF6khQZEamIE93cdANMU1fHeAMAAEzmrtu2qrWpWQ0NDaFuBUHAUBoAAHjt0zfeqcce+Iy6urpC3QoAP922oUxHm1tVX18f6lZwnQ0MDKi9cq/WpSzX691m/x1fv2qN6mp3B1ynqalJi4fmKCN+kTr6ze6WXpm2VJ2dnQHXOXDggG655RZjA9pR5eXlRmpalqXNpWXqGxgw0NUVaalpOnqkLeA6XV1dSktIkiOvUDsbdgXe2FU2pmXw5hUwjX106z363VPP6siRI6FuBQAATGNPfOM7uveee7jhdAZiKA0AALxGjDcwM3yfX/BmpZdffll/viJPkpSWMF/d/ZeM1Y6MjJT7Yp+GhoYCqmNZlrLi05QVn6ZDl08a6m6EI68o4KHvpUuXNG/ePGVmZurgwYOGOhuxadMmIzulLcuS02bXhrXr1NzaYqCzK1KTUwK+Ma2yslL23ELZc4tk7Q38Roar2ZevNX6zAACziPEGAACeEOM9czGUBgAAPiHGGwh/qYtS9NAnHyDGe5b59bZtuj1/sySpOD1b1R1mzlgeVbBmg/bt2xdQDcuytDZxqTLjU9TSd9pQZyM2FwZ+JnJtba1sNpsiIiI0b9489fX1GepOiouLU09Pj9xud0B1Tp04qcR5iXKW2LTDes1QdyMcJTZVVwd2VrVlWSovtGlD5hrVH2o01NkI+4ocVTKUBqY1YrwBAIA3iPGemRhKAwAAn336xjv1xa98hRhvIIzd+d7bifGeRUaju9PmLpAk2dKzVdVpdijtzA98J/I+q1bzFauEqDj1DF821NmI1ekZ2r9nb0A1LMtSeXm5JMlms6m2ttZEa2PWrVun1tZWv59//PhxrViaLklyFNtVWRPYAPnNyp2lAX+Nd++q1YbV2YqOjlakIjQ8PGyoOykxZo6ON5rdHQ7APGK8AQCAN4jxnnkYSgMAAJ/FRkbra08+qQ9+8IOhbgVAAIjxnj2uju6WpA1pGWo47v/wcyLOjcUBDSwvXLig4aNXzpFeFJOo7iFzO5ElKWbArf7+fr+fX1VVpZKSEkmS0+lURUWFqdbGagbyOaysrJSjxCZJWrZ0qdraj5pqTZKUt26Danf5P4gfHh7W4IVexURFS5I25qxX45HDptqTJC3VXJ08aTb6HYB5xHgDAABPiPGeeRhKAwAAvxSUlCgrP58YbyCMEeM9e1wd3S1JMVHRGna75QowKvpqCxPn61ir/0PQ6upqrY5PHfs4Kz5Vh/rMDheL1+Wprq7Or+e63W5duHBBc+bMkSSVlJSoqqrKZHvasmVLQENpy7K0ubRs7ONFycm62NtrojVJUnR0tAb7+/0eIjU2Nip/9dqxjx15haqoD/wc7auVpOdwrjQQBojxBgAA3iDGe2ZhKA0AAPz28KOPEuMNhDlivGe+N0d3j8pbnKHDZ143+lpL5i/S6dP+nQVtWZbWzlsy9nFWfKqaL58y1ZqkwHZzt7W1KTs7e+zjhIQEdXd3m2pNkrRs2TIdPuz/zuE9DQ1ak5k19rG92KaqXWYjvPPX56qx0b+zoC3LkjO/aOxje26hqvbuNtWaJKl05XqG0kCYIMYbAAB4gxjvmYOhNAAA8FtMbCwx3sAMQIz3zPbm6O5RJek52tl+wOhr2XILVVlZ6ddzq3ZaWh6zcOzjJbFJ6hw4N8UzfGcPoD/LsuR0Oq+5lpOTo46ODhOtjUlMTFSvH7ubBwYGFOmOUGTklV/zHcU27azy7887GUexze+hr2VZKi0oGfs4JSlZJ876dxPDZLIWLlF9tdmzvgEEDzHeAADAE2K8Zw6G0gAAICDEeAPhjxjvme3N0d2jbOnZ2vX6IaOvVV7k38DS5XLp9domRSli7FpkRISiFalht7lBRUL8HJ3pOO7Xcy3L0ubN134enU6nduzYYaK1MXa7XTU1vkdaNzQ0aNPGwmuuFeZv1O4G/+LKJ1PuLPV7KH30cItSFiy85trS1Bt0rvuCidYkjXzf6PgFbrIBwgQx3gAAwBvEeM8MDKUBAEDAiPEGwh8x3jPTZNHdkpSasECvd5v9d3v9yizV1foex9zU1KTFQ3PGXc+IX6SOfrO7pVemLVVnZ6fPzztw4IBWrlx5zTWn02k8Krq8vNyvmpZlqdTuuOZabGys+gYGTLUmSUpZlKKjR9p8fl5XV5fSEpLGXXfkFWpnw67AG7vKxrQM3qwCwggx3gAAwBvEeIc/htIAACBgxHgDMwMx3jPPZNHdo9IS5qu7/5Kx14uMjJT7Yp+GhoZ8ep5lWcqKTxt3PSs+TYcunzTVniTJkVfk89D30qVLmjdv3rjrmZmZOnjwoKnWJEmbNm3ya6e0ZVly2uzjrm9Yu07NrS0mWhuTmpzi841olZWVsucWjrtuzy2SZfhcafvytZwrDYQZYrwBAIAnxHiHv+hQNwAAmB1cLpfuvOsuNbcE9qbodL6byj04HOoWpuTuv3bINDw0pHem5xurf3WM99atW43VxfTyyCOP6Hcv/s73Jw5P7zcY3/z3I1ARcTFG603V34qYRGOvc3WM9ze/+U1jdRE6v962TZ+eILp7VHF6tqo7mvTnWeOHhf4qWLNB+/btU0FBgdfPsSxL+YlLx13PjE9RRXezsd4kaXOhTT+seNann1W1tbWy2WzjrkdERGjevHnq6+tTfHy8kf7i4uLU09Mjt9utiIgIz094w6kTJ5U4b/y/B84Sm3ZYrykrc7WR/iTJUWJTdXW13v72t3v9HMuydEuhY9z1DZlr9JlDjcZ6kyT7ihw9ZL2gT3ziE0brAlN55ZVX9I8PPChN8ffWHen932lv9J82F30vSXFz4rxe63a7VTBvmbHXvjrGm+NEAAD+unTpkm6+6R3qOnVGUhB+N4+JMlovYpq/l2hCf1+fHrz7w8bq3XXbVj1172/U0NCgjRs3GquL64OhNADgurh06ZLcMTF6qbY2oDqphvoJhsHq6R03d+m3gX3uvfHwo4/qnaWletvb3qbk5OSgvx6uvz/94Y+qecX381NdzaeC0I05pv9+zH33JqP1rsff31F3vvd2PfXM06qvr/dpqIjpZyy6+453T7rGlp6t55pqjA6lnfkjO5F9+f7ZZ9WqXOPfUEiIilPP8GVjvUnS6vQM7d+z16fnWJal8vLyCR+z2Wyqra1VWVmZifYkSevWrVNra6tWr/ZukHz8+HGtWJo+4WOOYrse/epj+psP3GOsv3Jnqf7fH37n01B6965aPfLO8T1ER0crUhEaHh5WVJSZN/kSY+boeKPZ3eGAJ08++aSe+s+fa0X68knXdCUnGH3N57P/xmi999z5LqP1fPXRrffopo9/UHfccYdWrVoV0l4AAOHpwIEDskfdoIfuHLk50fTv5udtZn8+JU3z9xKnqye+8R29festqtm1SzExZm88QHBN5w1nAAAgzBDjDcwMxHjPDJ6iuyVpQ1qGGo63Gn1d58Zin6KTL1y4oOGjk58bvSgmUd1DfSZaGxMz4FZ/f7/X66uqqlRSUjLhY06nUxUVFaZaG6vpy+ewsrJSjpLxO7kladnSpWprP2qqNUlS3roNqt3l/c0yw8PDGrzQq5ioie+L35izXo1HDptqT5K0VHN18qTZ6HfAk8hI3mYLFDHeAADAk9SUFD348f9FjHcY4v+WAQCAUVfHeAMIT1fHeCN8/XrbNt0+RXS3JMVERWvY7ZbL7Tb2ugsT5+tYq/dD0Orqaq2OnzwLJSs+VYf6zA4Xi9flqa6uzqu1brdbFy5c0Jw5cyZ8vKSkRFVVVSbb0+bNm30aSluWpc2lk+/UXpScrIu9vSZakzSyu3mwv9/roVFjY6PyV6+d9HFHXqEq6n0/R3sqJek5nCsNhKGrY7wBAAAm8/733aHm2n1qaGgIdSvwAUNpAABg3MOPPqovfuUr6urqCnUrAPx053tv19HmVtXX14e6FfhhLLp77gKPa/MWZ+jwmdeNvv6S+Yt0+vRpr9ZalqW185ZM+nhWfKpaLps9AsCX3dxtbW3Kzs6e9PGEhAR1d3ebak2SlJ6erkOHDnm9vqG+QWsysyZ93F5sU9WuahOtjclfn6vGRu/OgrYsS878okkft+cWqnqfdzcJeKt05XqG0kCY+ujWe/S7p57VkSNEmgIAgMk98cgXde8H/oqUtzDCUBoAABhHjDcwMxDjHb68ie4eVZKeo53tB4y+vi23UJWVlV6trdppaXnMwkkfXxKbpI6ByeO9/WH3oT/LsuR0Oqdck5OTo46ODhOtjZk/f756vdjdPDAwoChFTBkb7Ci2aWeVd39ebzmKbV4PfS3LUmnBxPHnkpSSlKzjZ8zeeJC1cInqq72PGAcwvRDjDQAAPElNStYDd36IGO8wwlAaAAAEBTHeQPgjxjt8eRPdPapk2Rrtet37XbneKC/ybmDpcrn0em2TohQx6ZrIiAhFK1LDbnODiYT4OTrTcdyrtZZlafPmqT+XTqdTO3bsMNHaGLvdrpoaz5HWDQ0N2rSxcMo1hfkbtbvB7E7kcmep10Ppo4dblLJg8hsPJGlp6g06133BRGuSRr5vdPwCN9UAYYoYbwAA4I27brqFGO8wwlAaAAAEDTHeQPgjxjv8+BLdLUlp85LU2X3WaA/rV2aprna3x3VNTU1aPDTxWc1Xy4hfpI5+s7ulV6YtVWdnp8d1Bw4c0MqVK6dc43Q6jUdFl5eXe1XTsiyV2h1TromNjVVff7+p1iRJKYtSdPRIm8d1XV1dSktI8rjOkVeonQ27Am/sKhvTMnhzCghjxHgDAABvEOMdPhhKAwCAoCHGG5gZiPEOL75Ed49anLBA3f2XjPUQGRkpd2+/hoaGplxnWZay4tM81suKT9OhyydNtSdJcuQVeRz6Xrp0SfPmzfNYKzMzUwcPHjTVmiRp06ZNXu2UtixLTpvd47oN69arubXFRGtjUpNTdO7c1DcLVFZWyp479U5uSbLnFsna6/lGBl/Yl6/lXGkgzBHjDQAAPCHGO3wwlAYAAEFFjDcQ/ojxDi++RHePKk7PVnVHk9E+CrLWa9++fVOusSxLaxOXeqyVGZ+ilr7TplqTJG0u9BwxXltbK5vN5rFWRESEEhMT1dfXZ6o9xcXFqaenR263e8p1p06cVOK8RI/1nCU27bBeM9WeJMlRYlNVVdWUayzLUnmh589h7ups1R9qNNWaJMm+IkeVDKWBsEaMNwAA8AYx3uGBoTQAAAg6YryB8EeMd3jwNbp7lC09W1WdZofSznzPO5H37tyl+Yr1WCshKk49w5dNtSZJWp2eoX0Ne6ZcY1mWysvLvapXUlKi2tpaE62NWbdunVpaJt/dfPz4ca1Ymu5VLUexXZU11aZak+TdudK1Nbu0YXW2x1pRUVGKVISGh4dNtafEmDk6tr/ZWD0AoUGMNwAA8AYx3tMfQ2kAABB0xHgDMwMx3tOfP9HdkrQhLUMNx1uN9uLcWDzlwPLChQtytZ/3ut6imER1D5nbiSxJsYNS/xRnLVdVVamkpMSrWk6nUxUVFaZaG6tZWVk56eOVlZVylHjehSxJy5YuVVv7UVOtSZLy1m1Q7a7JB/HDw8Ma6r6kmKhor+ptzFmvxiOHTbUnSVoWkaATJ04YrQng+iPGGwAAeEKM9/THUBoAAFwXxHgD4Y8Y7+nPn+huSYqJitaw2y2Xh6hoXyxMnK9jrZMPQaurq7U6PtXrelnxqTrUZ/Zc6eJ1eaqrq5vwMbfbrQsXLmjOnDle1SopKfEYZe2rzZs3TznYtyxLm0vLvK63KDlZF3t7TbQmSYqOjtZgf/+kQ6LGxkblr17rdT1HXqEq6j2fo+2LkvScKQf7AMIDMd4AAMAbxHhPbwylAQDAdUOMNxD+iPGevgYGBtR56IjP0d2j8hZn6PCZ1432tGT+Ip0+PfFZ0JZlae28JV7XyopPVcvlU6ZakzT1bu62tjZlZ3uOnR6VkJCg7u5uU61JktLT03Xo0KFJH2+ob9CazCyv69mLbaraZTbCO399rhobJz4L2rIsOfOLvK5lzy1U9b6JbxLwV+nK9R4jxgGEB2K8AQCAN4jxnr4YSgMAgOuGGG9gZiDGe3p6+eWXdZN9i9/PL0nP0c72AwY7kmy5hZPuUq3aaWl5zEKvay2JTVLHwDlTrUkaGYJO1p9lWXI6nT7Vy8nJUUdHh4nWxsyfP1+9E+xuHhgYUJQiFBnp/a/1jmKbdlaZ3TXsKLZNOvS1LEulBd7Fn0tSSlKyjp8xe+NB1sIlqq82e9Y3gNAhxhsAAHhCjPf05d3BTm+wLEutrWbPGQMAzA59fWbPgET4Go3xvu+++2S320PdDnzEEBLSGzHe9z2gW265hZtMppEf/vCH+tlDX5Uq/RuKlixboy/84f/qbzb9hbGeyots2ma9qptvvvma6y6XS6/XNilq4Q1e14qMiFC0IjXsdikqwsz91Qnxc3Sm4/iEj1mWpQceeMCnek6nUzt27NDdd99toj1Jkt1uV01Njd7ylrdcc72hoUGbNhb6VKswf6O+/C2z0bflzlI98qUv6MMf/vC4x44eblHKAu9vPJCkpak36Fz3BS2c79+O/zeLjIiQjl/Q4OCgYmJivH7e5z//eb322mu6ePGikT4k6dZbb9Xy5cuN1cP00t7eLrfBIxAwsdEY77vvvnvczxYAAEjTwKi7brpFTz30D3rssceUkZER6nZmrMzMTJ9u5vZ6KG1Zlt4ZmSBl5fnV2Gx3zp5vtN5Bo9Wk1dVm/7GOsa0yWu/CP//aaL2V7/I+hs8bluGv73NGq0n/WrXHcMXp7ZDhr0dyl7lz9yQpKjnBaL2JAzGnn4sXL+qFP/wh1G1gmnjo0Uf1ztJSffzzn9fC5GRJUpLhn0WXfmt2V9S8//UOo/W6DP9bYPrzF5mVNuH17377caOvg/B1522366lnn9ba9Ext+jPfdpN6Mmz4Z69ppn+WDxr4+zswOKD/6Ps3LUlJ0yX5N5ROm5ekzu6zAfdytfUrs1T3k2+Nu97U1KTFQ96d1Xy1jPhF6ug/p5Xxi0y0J0lambZUnZ2dSk9Pv+b6gQMHtHLlSp9qOZ1Offvb3zY6lC4vL9err746bihtWZZK7Q6fasXGxqqvv99Yb5KUsihFR4+0jbve1dWltIQkn+s58gq1s2GX3rX5bYE394aNaRlqaGhQcXGxV+vPnTun+vp6Pfvssx7XDp+75FXNliOt+vef/lj/+JFPTLnO1Wx2p7jp380XGv7d0vTvbqb58rvg//xymyIiIoLYDUZ9dOs9euZjH1BJwnKtWjp9bvSY7P+f/XXx2y8Yrbfg0duN1jP578GpiERjtSTpvOF/+0y/L2T6a2ua6d99Tf95Bwx/L2cb/tn2WcM/295ltJq01nA9E7/LXC3Q93Lqj8fod4Z6Qfh74uEv6h0f/ytZP/2NYqI936Bq+md5yuEWo/VM/7z87w8+GnCN0sP/oZ07d3o9mPb69nJ2SAMAAFNiY2P19Sef1D+wwxIIa9//+nf0tx/7CDvop4HfV1XoptK3BFxnccICdfd7N2TzRmRkpNy9/RoaGrrmumVZyor3/Rf+rPg0Hbp80lR7kiRHXtG4+OlLly5p3rx5PtfKzMzUwYNmbyHetGmTampqxl23LEtOm++JIxvWrVdzq9k3R1KTU3Tu3LXR6pWVlbLn+raTW5LsuUWy9u421dpIzeVrfTpX+plnntHWrVs1Z84cz//Fx3v134sv/07vuumdRv9cwGz2g898RXd/+hPEeAMAgEmlLkzWg/f8vT7zfbNpUbiWL/NjzpQGAAAhUVBSog35+Xpm27ZQtwLAT6kpKXro/gf10EMPhbqVWe9//vCi7vyLdwdcpzg9W9UdTQY6uqIga7327dt3zTXLsrQ2canPtTLjU9TSZzYnZnPh+DORa2trZbPZfK4VERGhxMREo8eWxMXFqaenZ1ws8KkTJ5U4z/c75Z0lNu2wXjPVniTJUWJTVVXVNdcsy1J5oe+fw9zV2ao/1GiqNUmSfUWOKn0YSj/99NO69dZbjfbw0h9f0Y2b/T/zHcC1Vi5N11+96zZ9479+GOpWAADANHbXTbfoyLEO7Tl0INStQAylAQBACD306KP67le+onNdXaFuBYCf7nzv7Wpvb1d9fX2oW5m1BgYHdOz0Cd2QEnjUmC09W1WdZofSzvzxO5H37tyl+Yr1uVZCVJx6hi+bak2StDo9Q/saro1ttCxL5eXlftUrKSlRba3ZIyTWrVunlpYru5uPHz+uFUvTp3jG5BzFdlXWVJtqTdLIudLjBvs1u7Rhte/HJkVFRSlSERoeHjbVnhJj5ujY/mav1p47d07R0dF+7ZSfTEdnp1ZlrFR0tNcnqAHwwke33qPfV+3QkWP+HVsBAABmhyce/qL+9gsPanCIlLdQYygNAABChhhvYGZ48sknde+99xLjHSKmorslaUNahhqOmz26ybmx+JqB5YULF+RqP+93vUUxieoeMrcTWZJiB6X+q85arqqqUklJiV+1nE6nKioqTLUmSXI4HKqsrBz7uLKyUo4S33chS9KypUvV1n7UVGuSpLx1G1S768ogfnh4WEPdlxQT5d8QdmPOejUeOWyqPUnSsogEnThxwuO6Z555RrfddpvR1/7l/2zT7bearQlgBDHeAADAE2K8pw+G0gAAIKRGY7x//cpzoW4FgJ9SU1P1yCOPEOMdIqaiuyUpJipaw263XG5zb+4vTJyvY61XhqDV1dVaHZ/qd72s+FQd7jN7rnTxujzV1dVJktxuty5cuKA5c+b4VaukpGRclHWgtmzZcs1g37IslTvL/K63KDlZF3t7TbQmSYqOjtZgf//YUKixsVH5q9f6Xc+RV6jX6neZak+SVJKec81gfzJPP/203vOe9xh9baK7geAhxhsAAHiDGO/pgaE0AAAIuYcefVRf++kP1HXhfKhbAeCnu+66ixjvEDAZ3T0qJzVdB890GqsnSUvmL9Lp0yNnQVuWpbXzlvhdKys+Vc2XT5lqTdK1u7nb2tqUne177PSohIQEdXd3m2pNkpSenq5Dhw6NfdxQ36Ds1Vl+17MX21S1y2yEd/76XDU2jpwFbVmWnPlFftey5xaqat9uU61JkkpXrh8XMf5mo9HdCQkJxl6X6G4g+IjxBgAA3iDGO/QYSgMAgJCLjY3V9x/+V937uftD3QqAABDjff2ZjO4eNTw8rKp2s+dK2/OKxnapVu20tDxmod+1lsQmqWPgnKnWJI0MQUf7syxLTqczoHo5OTnq6DA7HJk/f756e3s1MDCgKEUoMtL/X+cdxTbtrPK8a9gXzhL72NDXsiyVFfgXfy5JKUnJOn7G7I0HWQuXqKF66rO+ie4Gwhcx3gAAwBNivEOPoTQAAJgWitfnK3/NOmK8gTBGjPf1ZzK6e1TL+ZOq6TzkeaEPygpLZFmWXC6XjtU2KUoRfteKjIhQtCI1bDBiPCF+js50HJc0MlDdsiWwqGWn06nt27ebaG2M3W5XTU2NGhoatGljYUC1CvLytave7E7kModzbCjdfrhVixb4f+OBJC1NvUHnui+YaE3SyPeN+/iFKW+aIbobCF/EeAMAAG8Q4x1aDKUBAMC08bmP3E+MNxDmiPG+foIR3X1xsE+J61foeMQlYzUlaf3KLNXvrlNTU5PShvw7q/lqGfGL1NFvdrf0yrSl6uzs1MGDB5WRkRFQLafT6dX5xb4oLy+XZVmyLEuldkdAteLi4jQwMGCosxEpi1LU3nZUXV1dSk1YEHA9R16hdjaYPVd6Y1qGGhoaJnyM6G4g/BHjDQAAvEGMd+gwlAYAANNGbEwMMd7ADECM9/URjOjuZxsr9b6tt2vphix1D142VjcyMlKui32qqKhQVnzgQ/Ss+DQdunzSQGdXOPKK9OqrrxoZSmZmZurgwYMGurpi06ZNqqmpGYkXt9kDrrdh3Xo1t7YY6OyKlIWL9PLLL8ueG9hObkmy5xbJ2mt2N7d9+dpJz5UmuhuYGYjxBgAAnhDjHToMpQEAwLRCjDcQ/ojxvj6CEd393OFduvnmm+VwOlVt+Fzpgqz1euGFF7Q2cWnAtTLjU9TSd9pAV1dsLrTp2Weflc1mC7hWRESEEhMT1dfXZ6CzEXFxcerp6dHJ4yeUOC8x4HrOEpt2WK8Z6OwKR4lNzzzzjMoLA/8c5q7OVv2hRgNdXWFfkTPpDnaiu4GZgRhvAADgDWK8Q4OhNAAAmHaI8QbCHzHewRWs6O7Y7CWaP3++nE6nqjrM7vR15hdp/569mq/YgGslRMWpZ9jcTm5JWp2eoQP79qu8vNxIvZKSEtXW1hqpNWrdunVaOD/waGxJchTbVVlTbaTWqHJnqWqqq7VhdXbAtaKiohSpCA0b3O2YGDNHJ06cGHed6G5gZiHGGwAAeIMY7+uPoTQAAJh2iPEGZgZivIMnmNHdkpSfn6+GU0eN1ndu3KS+XnOD5EUxieoeMjuYXpaUqvz8fCO1nE6nKioqjNQa5XA4tMDQUHrZ0qVqazf7Nc5bt0GXui8qJsrMIHZjznodPG12qLRixYpxg2miu4GZhxhvAADgCTHe1x9DaQAAMC2Nxni/8NofQ90KAD+Nxnh/4QtfCHUrM04wo7slKSYmRhFLF8jldhurvzBxgdIWLjJWLys+VYf7ThmrJ0n23AJjZ0GXlJSoqqrKSK1RZWVlutxncLCfnKzunm5j9aKjo5W9YpWxIZAjr1CW6R37Tue4CG+iu4GZZzTG+9+e+q9QtwIAAKaxu266RUdeb1dLp9kbdjExhtIAAGDasucW6OLl3lC3ASAAxcXF6urqCnUbM0qwo7tHFdiKdfjccWOvIUlrV65W95CZc5az4lPVfNnsUNqRt0mWZRmplZCQMGFUdCAyMjJ0+swZY/XsxTZV7aoxVk+SNq3LV+ORw0Zq2XMLVdN5yEitUVu2bLnma0x0NzBz2XMLdb77QqjbAAAA05w9t5AI7+uEoTQAAJiWunt79IUffEvvLv/zULcCwE8ul0sf/vCH9dBDD4W6lRkl2NHdo5xOp3Ye2W/0dex5RWrpO22k1pLYJLX2nzVSa5Q9t2DcLtpAHT1q9o77+fPnq7fXzA1b67JzjJ8r7cgvUtXeOiO1UpKStfdEm5Fao7Kzs7Vnz56xj4nuBmamwaFBffwrn9VH3veBULcCAACmseaONv3m1ZeUk7E61K3MCgylAQDAtHT/Nx7V4//7XzQnPj7UrQDw0+OPP66bb75ZGRkZoW5lRgl2dPcoh8NhfJdqWUGxWgxFbg9FutWTKA0NDxupJ0nz5ibo7Flzg+6hoSHj50rbi4pVU1drpFbPxR5V795lpNaoso3Fqty720gtl9ulcxEDOt9/yUg9SYqMjFRkZOTYWfdEdwMz07/+++P65Pv/VilJyaFuBQAATFMul0t//8WH9d9ffFwRERGhbmdWYCgNAACmnRde+6MiIyNVVlAc6lYA+Km1tVW/+MUvdP/994e6lRnlekV3S9LixYt1PMLcMFCS1q1ao/Z+M3HurZdPacPGPO1tNnvm8KpVq9TZ2WmkVnR0tPGd12WOUmO7m3fV7VbPJbPHZKQkJavjpJnY9+aONm3Iz1Vl+wEj9UYVFRWpoaGB6G5ghqpr2q+a/Xt01023hroVAAAwjT3+i5/oPW+5SStuWBbqVmYNhtIAAGBa6e7t0We//1U9/uAXQt0KAD+5XC596EMf0i9/+UvuNjbsekV3j1q6IUvdvReNvVZUZKTckobdroBrHeo9qXe84x2qqDcbP+10Oo2cK33s2DHl5eXp4EGzQ/OijQXaZWin9J79+1RcUqzm1hYj9UalLkzWOQPnuO7YXa13vOMdquow+zksLS2VZVlEdwMz0ODQoD76pU/rP//lW6FuBQAATGPNHW369SvP6+N3/nWoW5lVGEoDAIBphdhuIPwR2x081yu6e5TD6ZS1x8wAdFRGXLI6B84FXOfIUJfe//73Gzu/eNSWLVuMDKW3b9+u0tJSJSYmqq+vz0BnI+Li4nSx96LcbndAdYaGhhQVG62ysjLtsF4z1N0Ie26hqvfXB1yncl+d7rjjDu072xF4U1cpLS1VZWUl0d3ADERsNwAA8ITY7tBhKA0AAKYNYruB8Edsd/Bcz+juUc4gDKWz4tPUcvl0wHUSc5bqhhtu0Nm+HgNdXZGVlaX9+/cHXMeyLDmdTtlsNtXWmv0c5qzJUcuR1oBq7D/YqKJNm+R0Oo3FgY8qLygxcq5059mTSk9PV3T6IrkM7K4ftWDBAh0/fpzobmCGIbYbAAB4g9ju0GEoDQAApgViu4HwR2x3cF3v6G5Jys/PV92hwAe0V1szJ03NfYENpc+4L6mw1CZJylizWmfOmzmnelRcXJz6+/sDqtHY2Kg1a9bI4XCooqLCUGcjHMW2gAfJO6ydcjgcWrZsmdo62g11NiI3K0d1B/cFVKPnUq9Sli+RJBXainXw7OsmWhuTkZGht771rUZrEt0NhA6x3QAAwBvEdocWQ2kAADAtENsNhD9iu4Prekd3S1JMTIzcc2KM7lJNiIrT+aFLAdU42HNMTqdT0shu7or6GhOtjSkuLlZdnf+x4P39/Zo7d64iIiJUUlKiqqoqg91Jm0vLVLUrsKF0ZU2VHA6HJCllcaou9vaaaE2SFB0VraHhYblc/n/fVO7dfc3XeGdbo6n2xmouWrTIaE2iu4HQIbYbAAB4Qmx36DGUBgAAIff7558nthsIc8R2B1cwort7Bi5PGd09amNRoZqOBhYV/WZJ0XPVPeT/OcvNl09dM7A0ERV9tbKysoDOld69e7dstpGd3AkJCeru7jbVmiQpfekyHW5pDqjGyTOnx4ayDocj4CH3m+VlrVXjkcN+P99qqB37GjscDtV0HjLVmqSRc6UDufHgzYjuBkKH2G4AAOANYrtDj6E0AAAIqZ7ubv3Lww8T2w2EMWK7gy8Y0d3/70DVlNHdoxwOhyrqzA4ss+JT1RJAhHd3orR8+XJJ0rp167Sv1ezA0uFwqLKy0u/nV1RUjO1ClqScnBx1dHSYaG1M4vz56vVzd/PZrrNasfJKooHD4dDOKv//vBNx5Bepaq//Q9/dTftUUFAgSUpJSdGp6MDi1N9s/fr12rNnj7F6RHcDoUFsNwAA8Aax3dMDQ2kAABBSn77vPn31iSeI7QbCGLHdwReM6O7fHqqeMrp71MhOZHM7SiUpe85iv4fSAxEurbLnjd0AERkZqdikeRoaHjbWX0JCgs6ePev386urq8d2Sksjn8Pt27ebaG2MvahYNXW1fj33tSprbBeyJBUWFmpXveHd5huLVbXPv+8bl9ulobhIxcbGjl1bnpej8/2Bxb5fLTIyUpGRkRocHDRSj+huIDSI7QYAAJ4Q2z19MJQGAAAh8/vnn1dkVJTs5eWhbgWAn1rbjhDbHWTBiO6+0NerMwluj9HdkrR48WIdO+//ruaJLI1NUnu/f0Pf1sun5Cx1XnNtU0mx9jYfNNHamFWrVqmzs9Ov5547d06JiYljHzudzoB2Xk+kzFGqyhr/drBb1VXXDKXj4uI0MDxkqjVJUkpSstpPHPPruc0dbcovLLjmmtPpVGX7AQOdXVFUVKSGhoaA6xDdDYQGsd0AAMAbxHZPHwylAQBASIzGdn/5e98LdSsA/ORyufSR+z5KbHeQBSO6+7mmGm2ylXi9/oaM5eruvWjs9SMjIuSWNOx2+fzcQ70nrxmoSiMDy4p6sxHjTqfTr3Oljx07pqysrGuuZWZm6uBBs0Pzoo0F2uXnTuk9+/dpw4YN11zLy89Tc2uLidbGpC5M1rnuCz4/b8fu6gm/xlUdZj+HpaWlAZ0dPorobuD6I7YbAAB4g9ju6YWhNAAACImx2O45c0LdCgA/ffeHT+pdb38nsd1B9uNnf6Xb3vp2ozV/02ipsLDQ6/UOp0PWHv8GoJPJiEtW58A5n593ZKhLRUVF11yz2+0BnV88kS1btvg1sNy+ffu4gWpERIQSExPV19dnqj3FxcXpYu9Fud1un543NDSkqNhoRUVFXXPd6XRqh/Wasf4kyZ5bqOr99T4/r3Jf3bjPYW5urvadNXsud2lpqZEd7ER3A9cfsd0AAMATYrunH6+zpbKzs3XetjGYvcxofLtPM/9ittw6s+UQIHPBljDtfXffHeoWME28/Nxz+uVPf6pf/vSnoW7Fe4Z/doSr9Tn81MOI5tYWffZfP6/Lly/rH//pkVC3M6PFx8XpU9941Ov1wyfOe1zzx9Y9+nxurtc1nU6n/ueJ/9BNzj/z+jmeZMWnqeXyaWXELfLpeYk5SxUfH3/NtaSkJJ3t6zHWmyRlZWVp//79Pj/Psix97GMfG3fdZrOptrZWZWVlJtqTJOWsyVHLkVZlZa72+jn7DzaqaNOmcdedTqc+98hn9TcfuMdYf+UFJXrm1Zd8/r7pPHtSS5YsueZaVFSUotMXyeV2KTLCzP31CxYs0IkTJwKq0dp2RH/Y/ifFL04y0tNsM9N+d/vav3wp1C3MCrsP7tXX/vMHGhoe0j3/zPEhXpvGv8/EhroBXGsaf69Imvb9PTjN680Gj/zZHaFuAdPE47/4iXbUVSvr1s2hbsWI6frzMjMz0+u1Xg+lS0pKVF1drUOHDvnVFABgduvr69MLf/hDqNvANNB94YI+96lP6cc//rFiY6fr/05hMt/6xjdD3QKmAZfLpY988mP60pe+pJSUlFC3M+O53W6fd8N68nBm5ridqFPJz8/XPx3yfUA7lTVz0vR/T1frrVrr9XPOuC+psNQ24WMZa1brzPkuo7vm4uLi1N/fr7i4OK+f09jYqDVr1oy77nA4VFFRYXQo7Si2qbKm2qeh9A5rpxwOx7jry5YtU1tHu7HeJCk3K0ef/4FvPzd6LvUqZfmSCR8rtBXr4NnXtT5luYn2JEkrVqzQiRMndMMNN/j1/P959jd66KGHlJeXZ6wnhKfvfOc7iowkkDDYBocG9Q9f+oy++73vKjExMdTtAACmmSNHjujS7/eGug1MA80dbfo/z2/Tz372M3ZJB1Gmj+8teD2UlkYG0yUl3p87BgDAqIsXLzKUhiTpM5/8pH7yf/6PysvLQ90K/PC9x78b6hYwDXz3h0/q1ve9V5/85CdD3Qquk5iYGLnnxBjdpZoQFafzQ5d8es7BnmP6c+fEZ/eOnCtdo/e85SYT7UmSiouLVVdXN+EQdyL9/f2aO3fuhG96lJSU6PHHHzfWmyRtLi3TV771dX3wzvd7/ZzKmip98O/+esLHUhan6mJvr+YlJBjpLzoqWkPDw3K5XF4P6yr37p70TQ2n06mdTzxjdCjtdDpVWVmp97znPX49/6U/vaLf/f5lRUf79PYKZqCnn36aNzyvg3/998f14D89ortJ4QIATKC2tlZPMZSe9Vxul/7+sUf00u9f1ooVK0LdDq7CLZwAAOC6+f3zzyshKoqBNBDGWtuO6FfPPKX77ycuc7bZWFSopqOtRmsmRc9V95D35yw3Xz415cCycu9uU61JksrKynw6V3r37t2y2SbeyZ2QkKDu7m5TrUmS0pcu0+GWZp+ec/LMaS1aNHFkusPhUNWuahOtjcnLWqvGI4e9Xm811E76NXY4HKrpNJve5u/Z4ZLU3t6uzDVZDKSB66Suab92dRxiIA0AAKb0+M9/ovfd834G0tMQQ2kAAHBd9HR364sPP6zvfe97oW4FgJ9GY7t/+atfsRtsFnI4HKqoMzuwzIpPVUvfaa/XdydKy5dPvEt23bp12tdqdmDpcDhUWVnp9fqKioopd1Xn5OSoo6PDRGtjEufPV29vr1drz3ad1YqVGZM+7nA4tLPK+z+vNxz5RaraW+f1+t1N+1RQUDDhYykpKToV3W+osxHZ2dnas2ePX8/9xS9+oTvu4MxC4HoYHBrSR7/8Gf33z/871K0AAIBprLnzqJ7a+bLuu+++ULeCCTCUBgAA18Wn77tP//bEE5ozZ06oWwHgp+/+8End/N5blZEx+VALM9fITmTvh4veyJ6z2Ouh9ECES6vseZPeEBEZGanYpHkaGh421l9CQoK6urq8Xl9VVTXpTmlp5HO4fft2E62NsRcVq6au1qu1FZU7pzzvq7CwULvqDe8231iiqn3efd+43C4NxUUqNjZ20jXL83J0vt+32PepREZGKioqSoODgz4/96WXXtKNN95orBcAk/vXf/+OPvWZh5SSkhLqVgAAwDTlcrv09198WL/69TZupJ+mGEoDAICgI7YbCH/EdmPx4sV6/dwpozWXxiapvf+sV2tbL5+Ss3TygaokbSop1t7mgyZaG7Ny5Up1dnZ6tfbcuXNKTEyc9HGn0+l3VPRkyhylqqzxbgd7ZU31lEPpuLg4DQwPmWpNkpSStFDtJ455tba5o035hQVTrnE6napsP2CgsysKCwvV0NDg03Pa29u1atUqoruB64DYbgAA4A1iu6c/htIAACCoiO0Gwh+x3Ri1ZOUKdfdeNFYvMiJCbknDbpfHtYd6T045UJVGBpYV9WYjxr0dJB87dkxr1qyZck1mZqaamppMtSZJKtpYoF1e7pTes3+fNmzYMOWavPw8Nbe2mGhtTOrCZJ3rvuBx3Y7dUw/NpZGvR1WH2RsPSktLfb5ZgOhu4PogthsAAHiD2O7wwFAaAAAEFbHdQPgjthujHE6HrD3eDUC9lRGXrM6Bcx7XtQ11qaioaMo1drvdp/OLvbFlyxavBpbbt2/3OFCNiIhQYmKi+vr6TLWnuLg49VzskdvtnnLd0NCQomKjFRUVNeU6p9OpHdZrxvqTJHtuoar313tcV7Wv3uPnMDc3V/vOmj2Xu7S01KezwyWiu4HrhdhuAADgCbHd4YOhNAAACBpiu4HwR2w3ruZ0Oo0PpbPi09Ry2fO50vNylik+Pn7KNUlJSerqM7eTW5KysrK0f/9+j+ssy/I4UJUkm82m2lqzn8O12WvVcqR1yjX7DzaqaNMmj7WcTqfXceDeKi8oUeVez2dVd5w9oSVLlky5JioqSjHpi+TyYne9txYsWKATJ054vZ7obuD6ILYbAAB4g9ju8MFQGgAABAWx3UD4I7Ybb5afn6+6Q54HtL5YMydNzX1TD6XPuC+pqMzmVb2M7NU6c77LRGtj4uLi1N/fP+WaxsZGj/Hd0hsR4xUVplqTJDmKbR4HyTusnV4NzZctW6a2jnZTrUmScrNyVHdw35Rrei71KnXFUq/qFdpLdPDs6yZaG7NixQqvB9NEdwPBR2w3AADwBrHd4YWhNAAACApiu4HwR2w33iwmJkbuOTFGd6kmRMXp/NClKdcc7Dkmh8PhVT2Hw6GK+hoTrY0pLi7W7t2T7/Tt7+/X3Llzvbp5o7i4WFVVVSbb0+bSMlXWTF2zsqZKdrvdq3opi1N1sbfXRGuSpOioaA0ND8vlmvz7pnLvbp++xjvbGk21J8n7s8MloruB64HYbgAA4Amx3eGHoTQAADCO2G4g/BHbjclsLCpU09Gpo6J9lRQ9V91Dk5+z3Hz5lFe7fKU34qe9iIr2RVlZ2ZRnDu/evVs2m3c7uRMSEtTd3W2qNUlS+tJlam5tmXLNyTOntWjRIq/qORwOVe0yG+Gdl7VWjUcOT/q41VDr9dfY4XCopvOQqdYkjZwd7s250kR3A8FHbDcAAPAGsd3hh6E0AAAwithuIPwR242pOBwOVdSZHVhmxaeqdYoI7+5Eafny5V7VWrdunfa1mh1YOhyOKQeWFRUVXu/ylaScnBx1dHSYaG1M4vz56p1kd/PZrrNasdL7xAOHw6GdVZ4HtL5w5Bepel/9pI/vbtqngoICr2qlpKToVPTUceq+ys7O1p49ezyuI7obCC5iuwEAgDeaW1uI7Q5DDKUBAIBRxHYD4Y/YbkxlZCdyndGa2XMWT3qu9ECES6vseV7fIBEZGanYpHkaGh421l9CQoK6uiY/p7qqqsrrndLSyOdw+/btJlobYy8qVk1d7YSPVVR6d570qMLCQu2qN7zbfGPJpDvYXW6XhuIiFRsb63W95Xk5Ot8/dey7LyIjIxUVFaXBwcEp1xHdjam43e5QtxD2iO0GAACejN5IT2x3+CFvCgBwXcTGxqqtuVmff/DBgOpM5zGn68SFULcwpcGWk+OulWWs1zuyi429BrHdM19cXJz+9z9/2ufnuc+bGxwEw0R/PwIR0/o7o/Wm6i8z+Qb97aa/NPZao7HdlhcxtpidFi9erNfPnTJac2lsktr7z074WOvlU3KW/rlP9TaVFGtv80EV5mww0Z4kaeXKlers7FR6evq4x86dO6fExESvazmdTn3rW9/SBz7wAWP9lTlK9eprO/SW8i3jHqusqdYHP/TXXteKi4vTwPCQsd4kKSVpodpPHJvwseaONuUXFvhUz+l0qvLVA3r7mk0GuhtRWFiohoYGFRdP/P9GRHdjKrm5ufri17+i+VP8W9AXH2P0NVtP1xitt+vlif8dnkze4pW6M3/8vzn+Go3tfuxuEpcAAP5JS0vT9qP71fPyZUnmfzfvv2GB0Xpx0/y9RFPe97Z3yrahwFi9x//t+7rtjtuJ7Q5D/CYFALguYmNj9fKLL+r06cmjOXF9tbe36+nPftvYUHo0trumqspIPUxPz7/4gjo7O0PdBq7y9+/YamwoTWw3vLVk5Qp1917U/IR5RupFRkTILWnY7VJUxLWBXocundTNPuzylUYGlhUvvWZ0KO10OmVZlrZu3XrN9WPHjmnNmjU+1crMzFRTU5Ox3iSpaGOBvv69b0/42J79+7Rhg2+fi7z8PDW3tigrc7WB7kakLkzWue4LWjj/2jfzKupqfNrJLY18PX76X78zOpQuLS2VZVmTDqWJ7sZUPve5z6m5ufn67pb+yEeu32u9ycDAgD57x4eNDaVHY7tffPUVI/UAALPT8uXL9T9Vf9D58+dD3Qre8Kc//Un1TY3GhtLNrS3a9tvfyLIsI/VwfTGUBgBcNwsXLtTChQtD3QbeEBUVZbQesd2zw9y5c5WdnR3qNnA1g8NjYrvhLYfTIWtPrW5y/pmxmhlxyeocOKeMuEXXXG8b7FJRUZFPtex2u37yze9LdxprT1u2bNETTzwxbii9fft2nweqERERSkxMVF9fn+Lj4430FxcXp4sXe+R2u6+5qWRoaEhRsdE+/9x3Op3aYb1mdChtzy1U9f76cd83lXvr9NjH7/WpVm5urvafNXuTVGlpqf7rv/5Ln/jEJyZ8/KWXXtKnPvUpo6+JmSMiIsLnG1TCWX+/2XPdie0GAJiSlpamtLS0ULeBNzQ3N+tIq5lNSmOx3dxIH7Y4UxoAAASM2G4g/I3Gdt9///2hbgVhwOl0ytoz8fnF/sqKT1PL5fFvVszLWebz4DYpKUlnL/eYak2SlJWVpf3794+7blmWz0NpSbLZbKqtNfs5XJu9Vi1HWq+5tv9go4o2+b6b2Ol0qrKm2lRrkqTygonPle44e0JLlizxqVZUVJSi05PlcrtMtacFCxboxIkTEz5GdDcQPKOx3XfffXeoWwEAANMYsd3hj6E0AAAIyGhs9/e+x9lvQLgithu+ys/PV92h8QPaQKyZk6bmvmuH0mfcl1RUZvOrXkb2ap0532WitTFxcXHjdgc2Njb6tTvS6XSqoqLCVGuSJEexbdwgeYe106+h+bJly9TW0W6qNUlSblaO6g7uu+Zaz6Vepa5Y6le9QnuJDp593URrY1asWDHhYJrobiA4RmO7//vn/x3qVgAAwDQ2Gtt93333hboVBIChNAAACAix3UD4I7YbvoqJiZF7TozRXaoJUXE6P3TpmmsHe47J4XD4Vc/hcKiivsZEa2OKi4u1e/eVnb79/f2aO3euXzdzFBcXq6qqymR72lxapsqaa2tW1lTJbrf7VS9lcaou9vaaaE2SFB0VraHhYblcV75vKvfuDuhrvLOt0VR7kq6cHf5mL730km688UajrwWA2G4AAOAZsd0zB0NpAADgN2K7gfBHbDf8lV9YoKajrZ4X+iApeq66h/rGPm6+fMqvXb7SG/HTE0RFB6KsrEyVlZVjH+/evVslJSV+1UpISFB3d7ep1iRJ6UuXqbm15ZprJ8+c1qJFiyZ5xtTsdruqdpmN8M7LWqvGI4fHPrYaav3+GjscDtV0HjLVmqSRs8Ov/hpLRHcDwUJsNwAA8Aax3TMHQ2kAAOAXYruB8EdsNwLhdDpVUWd2YJkVn6rWqyK8uxOl5cuX+1Vr3bp12tdqdmDpcDiuGVhWVFT4PVCVpJycHHV0dJhobUzi/PnqfWN389mus1qx0v8EBKfTqZ1VlZ4X+sCRX6TqffVjH+9u2qeCggK/aqWkpOhUdL/nhT7Izs7Wnj17rrlGdDdgHrHdAADAG8R2zywMpQEAgF+I7QbCH7HdCMTITuQ6ozWz5yweO1d6IMKlVfY8v2+YiIyMVGzSPA0NDxvrLyEhQV1dV86prqqqks3m35nX0sjncPv27SZaG2MvKlZNXa0kqaLSv/OkRxUWFmpXveHd5htLxnawu9wuDcVFKjY21u96y/NydL7/kueFXoqMjFRUVJQGBwfHrhHdDZhHbDcAAPCE2O6Zh6E0AADw2UvNu4ntBsIcsd0I1OLFi/X6uVNGay6NTVJ7/1lJUuvlU3KW+j9QlaRNJcXa23zQRGtjVq5cqc7OTknSuXPnlJiY6Hetyc4vDkSZo1SVNSM72CtrqgMaSsfFxWlgeMhUa5KklKSFaj9xTJLU3NGm/MKCgOo5nU5Vth8w0NkVhYWFamhokER0NxAMxHYDAABvENs98zCUBgAAPukZuKxHrW3EdgNhjNhumLJk5Qp19140Vi8yIkJuScNulw5dOhnQQFV6I2K83mzE+Ogg+dixY1qzZk1AtTIzM9XU1GSosxFFGwvGdkrv2b9PGzZsCKheXn7euHOqA5W6MFnnui+ooq7GyNe4qsPsjQelpaVjNwsQ3Q2YRWw3AADwBrHdMxNDaQAA4JOHXvixnvzZT4jtBsIYsd0wxeF0yNpTa7RmRlyyOgfOqW2wS0VFRQHVstvt15xfbMKWLVtkWZa2b98e8EA1IiJCiYmJ6uvrM9TdyO7mixd7NDg4qKjYaEVFRQVUz+l0aof1mqHuRthzC1W9v16Ve+sC/hzm5uZq/9lOQ52NKC0tHTs7nOhuwCxiuwEAgCfEds9cDKUBAIDXXmrerfjiTGK7gTDWdv4ksd0wxul0Gh9KZ8WnqeXyac3LWab4+PiAaiUlJens5R5DnY3IysrS/v37ZVlWwANVSbLZbKqtNfs5XJu9Vr975WUVbdoUcC2n0zkWB25KecHIudIdZ09oyZIlAdWKiopSdHqyXG6Xoe6kBQsW6MSJE0R3A4YR2w0AALxBbPfMxVAaAAB4hdhuIPy53G594rf/Rmw3jMnPz1fdof1Ga66Zk6b9fcdVVGYzUi8je7VOnj1tpNaouLg47du3L+D4bumNiPGKCgNdXeEotmnbM08bGZovW7ZMbR3tBrq6IjcrRzWNDUpdsdRIvUJ7iRpPm90tvWLFCv3oRz8iuhswhNhuAADgDWK7ZzaG0gAAwCvEdgPh7wfVz+s9H/8bYrthTExMjM7198rlMrdLNSEqTucHeuVwOIzUs9vtsvbsNlJrVHFxsQYGBozc3FFcXDwWFW1KubNUVbuqZbfbjdSLjotRd4+5HefRUdEaHh6WzWbmxgOHw6HK9gNGao1yOp165plniO4GDCG2GwAAeEJs98xHBhUAAPCI2G4g/LWdP6mnz+xV5f0/DnUrmGGycrL15f/4vhYtSPZqfc35gx7XuCRjA8tNmzbpn+77R53sOuPV+qjaRR7XHDt2TJmZmYG2JklKSEhQa2urnnzySa/Wuy4NeLXucl+fFi3y/GfxxsaNG/XFr39Zq1au8rjWfcq74fXwsEv5+fmBtiZpZCj9pY/9b0XI8xtXc570btf8kSNHtHDhQqK7AQNGY7sfu5vEJQAAMDliu2c+frsCAABTGo3trm7cE+pWAPhpNLb7V3/8LXcbw7jHHntMv/zlL9V16ZJX61d//B0e1/xZerpWrlwZYGcjCgoKtPXv/krt7V5GUJ/2PLRctWqV7rnnngA7u+JHP/qRXnzxRWP1JOk/f/afxmo98MAD+vnPf66zFy54XpzgXc07PvRX2rJlS2CNvSElJUX3ffVzampq8rj2ohdfX2nkXOmf/exngbYGzHqjsd0vvvpKqFsBAADT2Ghst2VZoW4FQcRQGgAATInYbiD8EduNYFq+fLkefPDBULcxqcjISP31X/91qNuYks1mM7YzPBjS0tL0yU9+MtRtTOnOO+8MdQsAJkBsNwAA8ITY7tmDoTQAAJgUsd1A+CO2GwAAhAKx3QAAwBvEds8eDKUBAMCEiO0Gwh+x3QAAIBQGXcPEdgMAAI+I7Z5dIkPdAAAAmJ6I7QbCH7HdAAAgFL7y6q+I7QYAAFMitnv2Yac0AAAYh9huIPwR2w0AAEKh4cQRNcR162t33x3qVgAAwDRGbPfsw1AaAABcg9huIPwR2w0AAEJh0DWs//Xij/TS7p2hbgUAAExjxHbPTsR3AwCAaxDbDYQ/YrsBAEAofOXVX+nBL3+e2G4AADApYrtnL3ZKAwCAMcR2A+GP2G4AABAKxHYDAABvENs9ezGUBgAAkojtBmYCYrsBAEAoENsNAAC8QWz37EZ8NwAAkERsNzATENsNAABCgdhuAADgCbHdYKc0AAAgthuYAYjtBgAAoUBsNwAA8Aax3WAoDQDALHehr5fYbiDMudwuYrsBAMB1NzA8RGw3AADw6HBHG7HdYCgNAMBs98jv/oPYbiDM/aD6BWK7AQDAdffV7duI7QYAAFNyuVz6/x57RL96mtju2Y4zpQEAmKXi4+O1s/Og4jatIrYbCGM9EYN6+sxe3X///aFuBQAAzBIRERE62ntWDXHdupvYbgAAMImFCxfq33/zc73vnvcT2w1FuN1ud6ibAAAAoTE8PCyXy6WYmJhQtwLAT263WwMDA4qLiwt1KwAAYBZxu90aHBxUbGxsqFsBAADT2MDAgGJiYtglDYbSAAAAAAAAAAAAAIDgIb4bAAAAAAAAAAAAABA0DKUBAAAAAAAAAAAAAEHDUBoAAAAAAAAAAAAAEDQMpQEAAAAAAAAAAAAAQcNQGgAAAAAAAAAAAAAQNAylAQAAAAAAAAAAAABBw1AaAAAAAAAAAAAAABA0DKUBAAAAAAAAAAAAAEHDUBoAAAAAAAAAAAAAEDQMpQEAAAAAAAAAAAAAQcNQGgAAAAAAAAAAAAAQNAylAQAAAAAAAAAAAABBw1AaAAAAAAAAAAAAABA0DKUBAAAAAAAAAAAAAEHDUBoAAAAAAAAAAAAAEDQMpQEAAAAAAAAAAAAAQcNQGgAAAAAAAAAAAAAQNAylAQAAAAAAAAAAAABBw1AaAAAAAAAAAAAAABA0DKUBAAAAAAAAAAAAAEHDUBoAAAAAAAAAAAAAEDQMpQEAAAAAAAAAAAAAQcNQGgAAAAAAAAAAAAAQNAylAQAAAAAAAAAAAABBw1AaAAAAAAAAAAAAABA0DKUBAAAAAAAAAAAAAEHDUBoAAAAAAAAAAAAAEDQMpQEAAAAAAAAAAAAAQcNQGgAAAAAAAAAAAAAQNAylAQAAAAAAAAAAAABBw1AaAAAAAAAAAAAAABA0DKUBAAAAAAAAAAAAAEHDUBoAAAAAAAAAAAAAEDQMpQEAAAAAAAAAAAAAQcNQGgAAAAAAAAAAAAAQNAylAQAAAAAAAAAAAABBw1AaAAAAAAAAAAAAABA0DKUBAAAAAAAAAAAAAEHz/wNjAU7ojHXjxAAAAABJRU5ErkJggg==\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 17
        },
        "id": "PxZS0WpKshib",
        "outputId": "f3609058-e9c6-4f6d-dca8-8d117f03ffe9"
      },
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<IPython.core.display.Javascript object>"
            ],
            "application/javascript": [
              "\n",
              "    async function download(id, filename, size) {\n",
              "      if (!google.colab.kernel.accessAllowed) {\n",
              "        return;\n",
              "      }\n",
              "      const div = document.createElement('div');\n",
              "      const label = document.createElement('label');\n",
              "      label.textContent = `Downloading \"${filename}\": `;\n",
              "      div.appendChild(label);\n",
              "      const progress = document.createElement('progress');\n",
              "      progress.max = size;\n",
              "      div.appendChild(progress);\n",
              "      document.body.appendChild(div);\n",
              "\n",
              "      const buffers = [];\n",
              "      let downloaded = 0;\n",
              "\n",
              "      const channel = await google.colab.kernel.comms.open(id);\n",
              "      // Send a message to notify the kernel that we're ready.\n",
              "      channel.send({})\n",
              "\n",
              "      for await (const message of channel.messages) {\n",
              "        // Send a message to notify the kernel that we're ready.\n",
              "        channel.send({})\n",
              "        if (message.buffers) {\n",
              "          for (const buffer of message.buffers) {\n",
              "            buffers.push(buffer);\n",
              "            downloaded += buffer.byteLength;\n",
              "            progress.value = downloaded;\n",
              "          }\n",
              "        }\n",
              "      }\n",
              "      const blob = new Blob(buffers, {type: 'application/binary'});\n",
              "      const a = document.createElement('a');\n",
              "      a.href = window.URL.createObjectURL(blob);\n",
              "      a.download = filename;\n",
              "      div.appendChild(a);\n",
              "      a.click();\n",
              "      div.remove();\n",
              "    }\n",
              "  "
            ]
          },
          "metadata": {}
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<IPython.core.display.Javascript object>"
            ],
            "application/javascript": [
              "download(\"download_347d0a21-14b5-4e1d-8ba1-07578cdbf254\", \"output.zip\", 147368)"
            ]
          },
          "metadata": {}
        }
      ],
      "source": [
        "#@title Download output files\n",
        "\n",
        "!zip -r output.zip output # first, zip files\n",
        "files.download(\"output.zip\")# Then download zipped file"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Optional: Creating a composite image\n",
        "If running SSDraw on multiple pdbs, you can use the code block below to combine your images into a composite image."
      ],
      "metadata": {
        "id": "yUC0X28sDRBJ"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "#@title Create composite image\n",
        "# create list and order all outputs by creation time.\n",
        "# then combine images together\n",
        "from PIL import Image\n",
        "\n",
        "img_names = [\"output/\"+img for img in os.listdir(\"output/\") if img != \"composite.png\"]\n",
        "img_names.sort(key = lambda x: os.path.getctime(x))\n",
        "\n",
        "imgs = [Image.open(img_name) for img_name in img_names]\n",
        "\n",
        "# combine SSDraw images into a composite image\n",
        "def combine_images(imgs):\n",
        "\n",
        "    height = min([img.height for img in imgs])\n",
        "    width = min([img.width for img in imgs])\n",
        "    pic = Image.new('RGB', (width, height*len(imgs)))\n",
        "\n",
        "    height_i = 0\n",
        "\n",
        "    for img in imgs:\n",
        "        x = img.crop((0,0,width,height))\n",
        "        pic.paste(x, (0, height_i))\n",
        "        height_i+=height\n",
        "    return pic\n",
        "\n",
        "combine_images(imgs).save(\"output/composite.png\")\n",
        "im = Image.open(\"output/composite.png\")\n",
        "im"
      ],
      "metadata": {
        "id": "PEfc9neUHZYN",
        "cellView": "form"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "vToOLIiTxHmp"
      },
      "source": [
        "# Instructions <a name=\"Instructions\"></a>\n",
        "**Quick start**\n",
        "1. Upload a pdb file and a sequence file in fasta format. To do this, click the folder icon on the left tab, and select the \"upload files\" option.\n",
        "2. Run the first cell block titled \"SSDraw\". This will download and install the necessary dependencies. It should take around 1-2 minutes.\n",
        "3. Fill out the arguments in the required section of \"Define Variables.\" Input the sequence file name, name of the protein used in the fasta file, pdb file name, and output file name into their respective prompts. Leave the other options unchanged and run this cell block.\n",
        "4. Run the two cell blocks \"Choose subregion to run SSDraw on\" and \"run SSDraw\". You should see an output image in the window.\n",
        "5. To download the results, run \"Download output files\".\n",
        "**Running on multiple PDBs**\n",
        "\n",
        "Run the first cell titled \"SSDraw\". Then fill in the required inputs for your first pdb, and run the next three cells: \"Define variables\", \"Choose subregion...\", and \"Run SSDraw\". That will generate an image for your first PDB. To run again for the next pdb, edit the input section and rerun the same three cells \"Define variables, \"Choose subregion...\", and \"Run SSDraw\". You do not need to rerun the first cell, since it is just for installing and importing the required libraries. Repeat this step for all PDBs.\n",
        "Each time you run SSDraw, the resulting image will be saved to the directory named \"output\". Run the final cell titled \"Download output files\" to download a zip file of the output directory.\n",
        "\n",
        "*Note*: If you are running SSDraw for multiple pdbs on a subregion of the alignment, you can skip rerunning the \"Choose subregion...\" cell. This will fix the subregion for all of your runs.\n",
        "\n",
        "**File formats:**\n",
        "\n",
        "The alignment file you use must be in fasta format.\n",
        "\n",
        "If you use a custom file to score residues, write the file in the following format: every line includes first the one-letter amino acid code of the residue, followed by a space, and then the score for that residue.\n",
        "\n",
        "**DSSP files:**\n",
        "\n",
        "Normally, SSDraw will generate a DSSP annotation from the PDB file, but if you have a DSSP file you would like to use, you can upload it and input the file name in Options.\n",
        "\n",
        "**Coloring:**\n",
        "\n",
        "SSDraw provides both solid and gradient coloring. If you choose solid, you can pick which solid color to use by either typing a color name contained in matplotlib's color library or a hexcode (with a preceding #). If you choose gradient, SSDraw will use the gradient color map selected in the options section. The image will be colored in based on a per-residue score (see Scoring). SSDraw also provides the option to color based on the mview color map, which colors residues by their chemical property.\n",
        "\n",
        "**Scoring**:\n",
        "\n",
        "There are two scoring methods built into SSDraw: scoring by conservation score or scoring by b-factor. There are also two ways to use a custom scoring method. First, you can upload a scoring file following this file format: every line includes first a one-letter amino acid code, followed by a space, and then the score of that residue. Another way to use a custom scoring is to edit the b-factors of the C-alpha carbons in your pdb file, then choose score by b-factor."
      ]
    }
  ],
  "metadata": {
    "colab": {
      "provenance": [],
      "authorship_tag": "ABX9TyMgNx6Sh2sWLShqPXpAI0yk",
      "include_colab_link": true
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}