HIV_analyze.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <font color = red> HIV Analysis\n",
    "</font> Here we combine HIV sequence data from [the Los Alamos National Laboratory HIV Sequence Database](https://www.hiv.lanl.gov/content/index) and immunological data to investigate HIV evolution across 13 individuals. This data is contained in the data/HIV/ directory.\n",
    "\n",
    "And we use some processed data from [here](https://github.com/bartonlab/paper-MPL-inference).\n",
    "- epitopes.csv $\\to$ epitopes.csv: information about epitopes \n",
    "- src-MPL-HIV.tar.gz/\\*-poly-seq2state.dat $\\to$ input/sequence/\\*-poly-seq2state.dat: processed sequences readable by MPL\n",
    "- src-MPL-HIV.tar.gz/Zanini-extended $\\to$ input/Zanini-extended: mutation matrix for HIV data\n",
    "- processed/\\*-index.csv $\\to$ notrait/processed/\\*-index.csv: information about all sites\n",
    "- interim/\\*-poly.csv $\\to$ notrait/interim/\\*-poly.csv: information about polymorphic sites, which will be modified to interim/*-poly.csv when finding escape sites"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contents\n",
    "- ### [Libraries and variables](#lib)\n",
    "- ### Data analyze\n",
    "    - #### [Find escape sites](#escape)\n",
    "    - #### [Write shell script](#mpl)\n",
    "    - #### [Deal with data](#deal)\n",
    "    - #### [Calculate Δs_ij](#sij)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a id='lib'></a> Libraries and variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "This notebook was prepared using:\n",
      "python version 3.10.9 | packaged by conda-forge | (main, Feb  2 2023, 20:26:08) [Clang 14.0.6 ]\n",
      "numpy version 1.24.2\n",
      "pandas version 1.5.3\n"
     ]
    }
   ],
   "source": [
    "print('This notebook was prepared using:')\n",
    "\n",
    "import os\n",
    "import sys\n",
    "print('python version %s' % sys.version)\n",
    "\n",
    "import numpy as np\n",
    "print('numpy version %s' % np.__version__)\n",
    "\n",
    "import pandas as pd\n",
    "print('pandas version %s' % pd.__version__)\n",
    "\n",
    "import re\n",
    "from math import isnan\n",
    "\n",
    "import HIV\n",
    "\n",
    "# GLOBAL VARIABLES\n",
    "\n",
    "NUC = ['-', 'A', 'C', 'G', 'T']\n",
    "\n",
    "ppt = ['700010040', '700010058', '700010077', '700010470', '700010607', '703010131', '703010159', \n",
    "       '703010256', '704010042', '705010162', '705010185', '705010198', '706010164']\n",
    "\n",
    "tags = ['700010040-3', '700010040-5', '700010058-3', '700010058-5', '700010077-3', '700010077-5',\n",
    "        '700010470-3', '700010470-5', '700010607-3', '700010607-5', '703010131-3', '703010131-5', \n",
    "        '703010159-3', '703010159-5', '703010256-3', '703010256-5', '704010042-3', '704010042-5', \n",
    "        '705010162-3', '705010162-5', '705010185-3', '705010185-5', '705010198-3', '705010198-5', \n",
    "        '706010164-3', '706010164-5']  \n",
    "\n",
    "# GitHub di rectories\n",
    "HIV_DIR = 'data/HIV'\n",
    "MPL_DIR = 'src'\n",
    "SIM_DIR = 'data/simulation'\n",
    "FIG_DIR = 'figures'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data analyze"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  <a id='escape'></a>Find escape sites\n",
    "escape site = polymorphic site which is nonsynonymous in the reading frame of its epitope.\n",
    "\n",
    "Steps to find escape sites:\n",
    "1. find all polymorphic sites within one epitope\n",
    "2. check if these polymorphic sites are nonsynonymous in the reading frame of its epitope\n",
    "3. Count the sites that meet the above conditions. \n",
    "    - \\>= 2, sites $\\to$ escape sites (or trait sites), epitope $\\to$ escape group (or trait group);\n",
    "    - <  2, sites $\\to$ special sites, epitope $\\neq$ escape group (or trait group). \n",
    "\n",
    "Special sites: nonsynonymous mutation in the reading frame of their epitope but not escape sites because the number of escape sites in its epitope is less than 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "++ 700010040-3 has 3 (+1) escape groups, they are [18 19 20 21],[257 258 259],[286 287 288 290 291], the special sites are [245 246]\n",
      "== 700010040-5 has no escape group, the special sites are: [31 32],[39],[104]\n",
      "-- 700010058-3 has 1 escape groups, they are [60 61 62 63]\n",
      "== 700010058-5 has no escape group, the special sites are: [13],[16 17]\n",
      "-- 700010077-3 has 6 escape groups, they are [28 31 32 34 35],[82 83 84 85],[87 88 89],[122 123 124 125],[149 150 151 152 153 154 156 157 158 159 160 161 162 163 164 165 166 167\n",
      " 168 169 170 171],[186 187 188 189 190]\n",
      "== 700010077-5 has no escape group, the special sites are: [1],[4 5]\n",
      "++ 700010470-3 has 3 (+2) escape groups, they are [0 1 2 3],[302 303 304 305],[352 353 355], the special sites are [40],[334]\n",
      "++ 700010470-5 has 3 (+2) escape groups, they are [23 24 25],[78 79 81],[149 150 151 152], the special sites are [8 9],[124]\n",
      "-- 700010607-3 has 1 escape groups, they are [198 199 200 201 202]\n",
      "   700010607-5 has no triat site\n",
      "\tmutant at site 0 in codon for CH131-3 that does not terminate in alignment\n",
      "++ 703010131-3 has 1 (+5) escape groups, they are [618 619 620 621 624], the special sites are [45 46],[47],[76],[189],[560 561]\n",
      "++ 703010131-5 has 1 (+1) escape groups, they are [14 15 16 19], the special sites are [259]\n",
      "-- 703010159-3 has 2 escape groups, they are [370 371 372 374 375 376],[417 418 420 421 422 423 424]\n",
      "   703010159-5 has no triat site\n",
      "++ 703010256-3 has 4 (+2) escape groups, they are [28 29 30 31 32 33],[345 346 347 348 349],[362 363 364 366],[423 424 427], the special sites are [45 46],[316 317]\n",
      "++ 703010256-5 has 1 (+2) escape groups, they are [169 170 171 172 174 175 176 177], the special sites are [42],[208]\n",
      "++ 704010042-3 has 2 (+2) escape groups, they are [200 201 202 203 204 205 207 208 210],[660 661 662 663 664 665], the special sites are [95 96],[765 766]\n",
      "== 704010042-5 has no escape group, the special sites are: [29 30],[125]\n",
      "++ 705010162-3 has 2 (+2) escape groups, they are [290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307\n",
      " 308 309 310 311 312 313 314 315 316 317],[334 335 337], the special sites are [33 34],[327 328]\n",
      "++ 705010162-5 has 3 (+1) escape groups, they are [6 7 8],[30 31 33],[38 39 40 42 43], the special sites are [155]\n",
      "   705010185-3 has no triat site\n",
      "-- 705010185-5 has 1 escape groups, they are [2 3 4 5 6]\n",
      "-- 705010198-3 has 1 escape groups, they are [194 195 196]\n",
      "   705010198-5 has no triat site\n",
      "++ 706010164-3 has 3 (+1) escape groups, they are [86 87 88 89],[182 183 184 185 186 187 188],[453 454 455 456 457 458], the special sites are [463]\n",
      "== 706010164-5 has no escape group, the special sites are: [44],[50 51],[170]\n"
     ]
    }
   ],
   "source": [
    "import importlib\n",
    "importlib.reload(HIV)\n",
    "    \n",
    "min_n = 2\n",
    "\n",
    "for tag in tags:\n",
    "    HIV.find_trait_site(tag,min_n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "tag|escape group|escape sites|special sites\n",
    ":----:|:----:|:---|:-----\n",
    "700010040-3|3|[18,19,20,21]<br>[257,258,259]<br>[286,287,288,290,291]|[245,246]\n",
    "700010040-5|0|                                                       |[31,32], [39], [104]\n",
    "700010058-3|1|[60,61,62,63]                                          |\n",
    "700010058-5|0|                                                       |[13], [16,17]\n",
    "700010077-3|6|[28,31,32,34,35]<br>[82,83,84,85]<br>[87,88,89]<br>[122,123,124,125]<br>[149,150,151,152,153,154,156,157,158,159,160,<br>161,162,163,164,165,166,167,168,169,170,171]<br>[186,187,188,189,190]|\n",
    "700010077-5|0|                                                       |[1], [4,5]\n",
    "700010470-3|3|[0,1,2,3]<br>[302,303,304,305]<br>[352,353,355]        |[40], [334]\n",
    "700010470-5|3|[23,24,25]<br>[78,79,81]<br>[149,150,151,152]          |[8,9], [124]\n",
    "700010607-3|1|[198,199,200,201,202]                                  |[]\n",
    "700010607-5|0|no mutation related to epitope                         |\n",
    "703010131-3|1|[618,619,620,621,624]                                  |[[45,46], [47], [76], [189], [560,561]\n",
    "703010131-5|1|[14,15,16,19]                                          |[259]\n",
    "703010159-3|2|[370,371,372,374,375,376]<br>[417,418,420,421,422,423,424]|\n",
    "703010159-5|0|no mutation related to epitope                         |\n",
    "703010256-3|4|[28,29,30,31,32,33]<br>[345,346,347,348,349]<br>[362,363,364,366]<br>[423,424,427]|[45,46], [316,317]\n",
    "703010256-5|1|[169,170,171,172,174,175,176,177]                      |[42], [208]\n",
    "704010042-3|2|[200,201,202,203,204,205,207,208,210]<br>[660,661,662,663,664,665]|[95,96], [765,766]\n",
    "704010042-5|0|                                                       |[29,30], [125]\n",
    "705010162-3|2|[6,7,8]<br>[30,31,33]<br>[38,39,40,42,43]              |[33,34], [327,328]\n",
    "705010162-5|3|[6,7,8]<br>[30,31,33]<br>[38,39,40,42,43]              |[155]\n",
    "705010185-3|0|no mutation related to epitope                         |\n",
    "705010185-5|1|[2,3,4,5,6]                                            |\n",
    "705010198-3|1|[194,195,196]                                          |\n",
    "705010198-5|0|no mutation related to epitope                         |\n",
    "706010164-3|3|[86,87,88,89]<br>[182,183,184,185,186,187,188]<br>[453,454,455,456,457,458]|[463]\n",
    "706010164-5|0|                                                       |[44], [50,51], [170]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check if the TF sequence and distance is correct "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "no problem\n"
     ]
    }
   ],
   "source": [
    "import importlib\n",
    "importlib.reload(HIV)\n",
    "\n",
    "TF_correct = True\n",
    "for tag in tags:\n",
    "    traitsite = HIV.read_file('traitsite/traitsite-'+tag+'.dat')\n",
    "    TFseq     = HIV.read_file('traitseq/traitseq-'+tag+'.dat')\n",
    "    traitdis  = HIV.read_file('traitdis/traitdis-'+tag+'.dat')\n",
    "    df = pd.read_csv('%s/interim/%s-poly.csv' %(HIV_DIR,tag), comment='#', memory_map=True)\n",
    "    for i in range(len(traitsite)):\n",
    "        for j in range(len(traitsite[i])):\n",
    "            n_ij   = df[df['polymorphic_index'] == traitsite[i][j]]\n",
    "            TF     = n_ij.iloc[0].TF\n",
    "            TF_nuc = NUC.index(TF)\n",
    "            if TFseq[i][j] != TF_nuc:\n",
    "                TF_correct = False\n",
    "                print('there is something wrong with TF sequence for %s'%tag)\n",
    "        \n",
    "        if len(traitsite[i]) != len(traitdis[i])+1:\n",
    "            TF_correct = False\n",
    "            print('there is something wrong with the trait distance for %s'%tag)\n",
    "                \n",
    "if TF_correct:\n",
    "    print('no problem')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check if there are some synonymous mutations on trait sites"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CH077-3 35 C\n",
      "CH077-3 161 A\n",
      "CH162-3 291 T\n",
      "CH162-3 295 T\n"
     ]
    }
   ],
   "source": [
    "for tag in tags:\n",
    "    traitsite = HIV.read_file('traitsite/traitsite-'+tag+'.dat')\n",
    "    df = pd.read_csv('%s/interim/%s-poly.csv' %(HIV_DIR,tag), comment='#', memory_map=True)\n",
    "    for i in range(len(traitsite)):\n",
    "        for j in range(len(traitsite[i])):\n",
    "            n_ij   = df[df['polymorphic_index'] == traitsite[i][j]]\n",
    "            for a in range(len(n_ij)):\n",
    "                if n_ij.iloc[a].escape == False and n_ij.iloc[a].nucleotide != n_ij.iloc[a].TF:\n",
    "                    print(f'CH{tag[-5:]} {n_ij.iloc[a].polymorphic_index} {n_ij.iloc[a].nucleotide}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a new input file to store allele information that does not contribute to trait term\\\n",
    "Previous : only wild type allele\\\n",
    "Special case: some synonymous mutant allele and wild type allele"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "special_tags = ['700010077-3', '705010162-3']\n",
    "\n",
    "# case 1: CH077-3 35C and 161A\n",
    "tag = '700010077-3'\n",
    "traitseq  = HIV.read_file('traitseq/traitseq-'+tag+'.dat')\n",
    "traitsite = HIV.read_file('traitsite/traitsite-'+tag+'.dat')\n",
    "f = open(\"%s/input/traitseq/traitallele-%s.dat\"%(HIV_DIR,tag),'w')\n",
    "for i in range(len(traitsite)):\n",
    "    line_content = [\n",
    "        f'{traitseq[i][j]}/{NUC.index(\"C\")}' if traitsite[i][j] == 35 else\n",
    "        f'{traitseq[i][j]}/{NUC.index(\"A\")}' if traitsite[i][j] == 161 else\n",
    "        str(traitseq[i][j])\n",
    "        for j in range(len(traitsite[i]))\n",
    "    ]\n",
    "    \n",
    "    line = '\\t'.join(line_content)\n",
    "    f.write(line + '\\n')\n",
    "f.close()\n",
    "\n",
    "# case 2: CH162-3 291T and 295T\n",
    "tag = '705010162-3'\n",
    "traitseq  = HIV.read_file('traitseq/traitseq-'+tag+'.dat')\n",
    "traitsite = HIV.read_file('traitsite/traitsite-'+tag+'.dat')\n",
    "f = open(\"%s/input/traitseq/traitallele-%s.dat\"%(HIV_DIR,tag),'w')\n",
    "for i in range(len(traitsite)):\n",
    "    line_content = [\n",
    "        f'{traitseq[i][j]}/{NUC.index(\"T\")}' if traitsite[i][j] == 291 else\n",
    "        f'{traitseq[i][j]}/{NUC.index(\"T\")}' if traitsite[i][j] == 295 else\n",
    "        str(traitseq[i][j])\n",
    "        for j in range(len(traitsite[i]))\n",
    "    ]\n",
    "    \n",
    "    line = '\\t'.join(line_content)\n",
    "    f.write(line + '\\n')\n",
    "f.close()   \n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id='mpl'></a> Write shell script\n",
    "\"src/HIV.sh\" : running MPL code for all HIV data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = open(\"src/HIV.sh\",'w')\n",
    "f.write('g++ main.cpp inf.cpp io.cpp -O3 -mcpu=apple-a14 -std=c++11 -lgsl -lgslcblas -o mpl\\n')\n",
    "for tag in tags:\n",
    "    if tag not in special_tags:\n",
    "        f.write('./mpl -d ../data/HIV -i input/sequence/%s-poly-seq2state.dat '%tag)\n",
    "        f.write('-o output/sc-%s.dat -g 10 -m input/Zanini-extended.dat -rr 1.4e-5 '%tag)\n",
    "        f.write('-e input/traitsite/traitsite-%s.dat -es input/traitseq/traitseq-%s.dat '%(tag,tag))\n",
    "        f.write('-ed input/traitdis/traitdis-%s.dat\\n'%tag)\n",
    "\n",
    "f.write('g++ main.cpp inf_special.cpp io.cpp -O3 -mcpu=apple-a14 -std=c++11 -lgsl -lgslcblas -o mpl\\n')\n",
    "# special case (escape sites contain nonsysnmous mutations)\n",
    "for tag in special_tags:\n",
    "    f.write('./mpl -d ../data/HIV -i input/sequence/%s-poly-seq2state.dat '%tag)\n",
    "    f.write('-o output/sc-%s.dat -g 10 -m input/Zanini-extended.dat -rr 1.4e-5 '%tag)\n",
    "    f.write('-e input/traitsite/traitsite-%s.dat -es input/traitseq/traitallele-%s.dat '%(tag,tag))\n",
    "    f.write('-ed input/traitdis/traitdis-%s.dat\\n'%tag)\n",
    "\n",
    "f.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"src/HIV_notrait.sh\" : running previous MPL code (no trait group) for all HIV data\\\n",
    "In this case, we use the same regularization value (Here gamma = g/N). \\\n",
    "No recombination part since $R_{fl} = 0$ for individual selection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = open(\"src/HIV_notrait.sh\",'w')\n",
    "f.write('g++ notrait/main.cpp notrait/inf.cpp notrait/io.cpp -O3 -mcpu=apple-a14 -lgsl -lgslcblas -o mpl\\n')\n",
    "for tag in tags:\n",
    "    f.write('./mpl -d ../data/HIV -i input/sequence/%s-poly-seq2state.dat '%tag)\n",
    "    f.write('-o notrait/output/sc-%s.dat -g 1e4 -N 1e3 -m input/Zanini-extended.dat\\n'%tag)\n",
    "f.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"src/HIV_noR.sh\" : running new MPL code but no recombiination term part for all HIV data\n",
    "We set recombination rate rr to 0 to remove the recombination flux $R_{fl}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = open(\"src/HIV_noR.sh\",'w')\n",
    "f.write('g++ main.cpp inf.cpp io.cpp -O3 -mcpu=apple-a14 -std=c++11 -lgsl -lgslcblas -o mpl\\n')\n",
    "for tag in tags:\n",
    "    if tag not in special_tags:\n",
    "        f.write('./mpl -d ../data/HIV -i input/sequence/%s-poly-seq2state.dat '%tag)\n",
    "        f.write('-o noR/output/sc-%s.dat -g 10 -m input/Zanini-extended.dat -rr 0 '%tag)\n",
    "        f.write('-e input/traitsite/traitsite-%s.dat -es input/traitseq/traitseq-%s.dat '%(tag,tag))\n",
    "        f.write('-ed input/traitdis/traitdis-%s.dat\\n'%tag)\n",
    "\n",
    "# special case (escape sites contain nonsysnmous mutations)\n",
    "f.write('g++ main.cpp inf_special.cpp io.cpp -O3 -mcpu=apple-a14 -std=c++11 -lgsl -lgslcblas -o mpl\\n')\n",
    "for tag in special_tags:\n",
    "    f.write('./mpl -d ../data/HIV -i input/sequence/%s-poly-seq2state.dat '%tag)\n",
    "    f.write('-o noR/output/sc-%s.dat -g 10 -m input/Zanini-extended.dat -rr 0 '%tag)\n",
    "    f.write('-e input/traitsite/traitsite-%s.dat -es input/traitseq/traitallele-%s.dat '%(tag,tag))\n",
    "    f.write('-ed input/traitdis/traitdis-%s.dat\\n'%tag)\n",
    "\n",
    "f.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id='deal'></a> Deal with data\n",
    "Dealing with mpl results (Normalize the selection coefficients such that the TF nucleotide at each site has a selection coefficient of zero) and record information in these 2 csv files.\n",
    "- /analysis/'tag'-analyze.csv \n",
    "    - coefficients in 2 cases (old: without escape part; new: with escape part)\n",
    "    - allele frequencies\n",
    "- /group/escape_group-'tag'.csv \n",
    "    - escape coefficients\n",
    "    - escape frequencies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Successfully create the analysis and escape group file\n"
     ]
    }
   ],
   "source": [
    "import importlib\n",
    "importlib.reload(HIV)\n",
    "\n",
    "for tag in tags:\n",
    "    HIV.analyze_result(tag,True)\n",
    "    HIV.analyze_result(tag,False)   # without recombination flux\n",
    "\n",
    "print('Successfully create the analysis and escape group file')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id='sij'></a> Calculate Δs_ij\n",
    "To quantify the linakge between different variants ($\\Delta \\hat s_{ij}$), we modify the input sequence by converting mutant variant $i$ into wild variant. \n",
    "\n",
    "$\\Delta \\hat s_{ij} = \\hat s_j - \\hat s_j^{\\setminus i}$\n",
    "- $\\hat s_j$ ---  Selection coefficient for variant j\n",
    "- $\\hat s_j^{\\setminus i}$ ---Selection coefficient for variant j calculated by using WT to replace variant i \n",
    "\n",
    "The modifications for individual loci and escape groups are different:\n",
    "- individual loci part\n",
    "   - replace all mutant variant $i$ by its transmitted/founder (TF) nucleotide for all sequences at all time points. \n",
    "- escape part\n",
    "   - change all mutant variants within one epitope into wild type."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "import importlib\n",
    "importlib.reload(HIV)\n",
    "\n",
    "# create the path\n",
    "for tag in tags:\n",
    "    # input sequence\n",
    "    path = 'data/HIV/input/sequence/'+ tag\n",
    "    folder = os.path.exists(path)\n",
    "    if not folder:\n",
    "        os.makedirs(path)\n",
    "    # output inferred coefficients\n",
    "    path = 'data/HIV/output/'+ tag\n",
    "    folder = os.path.exists(path)\n",
    "    if not folder:\n",
    "        os.makedirs(path)\n",
    "\n",
    "g = open(\"src/HIV_sij.sh\",'w')\n",
    "g.close()\n",
    "\n",
    "for tag in tags:\n",
    "    HIV.modify_seq(tag,special_tags) # Add content to shell script here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dealing with these data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "import importlib\n",
    "importlib.reload(HIV)\n",
    "\n",
    "for tag in tags:\n",
    "    HIV.cal_sij(tag)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}