0_create_gtf.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ee97c6cf-ecb4-4d43-9fe1-ee243243f193",
   "metadata": {},
   "source": [
    "# Generating a custom GTF file\n",
    "example HIV"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b9483296-dce9-4eec-8130-ca9631ab6ac3",
   "metadata": {},
   "source": [
    "### Data from paper (linked below) in computer-readable form:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b0c0654-32db-4b2a-9971-9b58f262c8cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "#manually entered donor and acceptor site positions, as well as isoform names <-> splice site information\n",
    "\n",
    "#https://retrovirology.biomedcentral.com/articles/10.1186/s12977-020-00533-1/figures/2\n",
    "donor_sites =  {\n",
    "        \"1\" : 743,\n",
    "        \"1c\" : 747,\n",
    "        \"1d\" : 793,\n",
    "        \"2\" : 4962,\n",
    "        \"2c\" : 4973,\n",
    "        \"2d\" : 5006,\n",
    "        \"2b\" : 5058,\n",
    "        \"3\" : 5463,\n",
    "        \"4\" : 6044,\n",
    "        \"4a\" : 6048,\n",
    "        \"4b\" : 6051\n",
    "    }\n",
    "\n",
    "acceptor_sites =  {\n",
    "        \"1\" : 4913,\n",
    "        \"1b\" : 4919,\n",
    "        \"1c\" : 4922,\n",
    "        \"2\" : 5390,\n",
    "        \"3\" : 5777,\n",
    "        \"4c\" : 5936,\n",
    "        \"4a\" : 5954,\n",
    "        \"4b\" : 5960,\n",
    "        \"5\" : 5976,\n",
    "        \"5a\" : 5980,\n",
    "        \"5b\" : 5983,\n",
    "        \"5c\" : 5992,\n",
    "        \"7c\" : 8345,\n",
    "        \"7\" : 8369\n",
    "    }\n",
    "\n",
    "D1_in_reference = 290\n",
    "offset = donor_sites[\"1\"] - D1_in_reference\n",
    "print(\"Detected offset of\", offset, \"between own reference and published data\")\n",
    "for site, value in donor_sites.items():\n",
    "    donor_sites[site] = value - offset\n",
    "    \n",
    "for site, value in acceptor_sites.items():\n",
    "    acceptor_sites[site] = value - offset\n",
    "\n",
    "isoforms =  {\n",
    "        \"unspliced1\" : [],\n",
    "        \"unspliced2\" : [],\n",
    "\n",
    "        \"Vif2\" : [\"D1A1\"],\n",
    "        \"Vif1\" : [\"D1A1\",\"D4A7\"],\n",
    "        \"Vif3\" : [\"D1cA1b\"],\n",
    "        \"Vif4\" : [\"D1cA1\"],\n",
    "        \"Vif5\" : [\"D1A1b\"],\n",
    "    \n",
    "        \"Vif6\" : [\"D1cA1b\", \"D4A7\"],\n",
    "        \"Vif7\" : [\"D1cA1c\"],\n",
    "        \"Vif8\" : [\"D1A1b\", \"D4A7\"],\n",
    "        \"Vif9\" : [\"D1cA1\", \"D4A7\"],\n",
    "\n",
    "        \"Vpr3\" : [\"D1A2\"],\n",
    "        \"Vpr4\" : [\"D1A1\", \"D2A2\"],\n",
    "        \"Vpr1\" : [\"D1A2\", \"D4A7\"],\n",
    "        \"Vpr2\" : [\"D1A1\", \"D2A2\", \"D4A7\"],\n",
    "\n",
    "        \"Env1\" : [\"D1A5\"],\n",
    "        \"Env9\" : [\"D1A2\", \"D3A5\"],\n",
    "        \"Env2\" : [\"D1A4b\"],\n",
    "        \"Env3\" : [\"D1A4a\"],\n",
    "        \"Env17\" : [\"D1A5a\"],\n",
    "        \"Env5\" : [\"D1A1\", \"D2A5\"],\n",
    "        \"Env13\" : [\"D1A1\", \"D2A2\", \"D3A5\"],\n",
    "        \"Env18\" : [\"D1cA4b\"],\n",
    "        \"Env10\" : [\"D1A2\", \"D3A4b\"],\n",
    "        \"Env19\" : [\"D1cA5\"],\n",
    "    \n",
    "        \"Env11\": [\"D1A2\", \"D3A4a\"],\n",
    "        \"Env14\" : [\"D1A1\", \"D2A2\", \"D3A4b\"],\n",
    "        \"Env15\" : [\"D1A1\", \"D2A2\", \"D3A4a\"],\n",
    "        \"Env12\" : [\"D1A2\", \"D3A4c\"],\n",
    "        \"Env6\" : [\"D1A1\", \"D2A4b\"],\n",
    "        \"Env7\" : [\"D1A1\", \"D2A4a\"],\n",
    "        \"Env4\" : [\"D1A4c\"],\n",
    "        \n",
    "        \"Env8\" : [\"D1A1\", \"D2A4c\"],\n",
    "        \"Env20\" : [\"D1A2\", \"D3A5a\"],\n",
    "\n",
    "\n",
    "        \"Tat1\" : [\"D1A3\", \"D4A7\"],\n",
    "        \"Tat2\" : [\"D1A1\", \"D2A3\", \"D4A7\"],\n",
    "        \"Tat5\" : [\"D1A3\"],\n",
    "        \"Tat3\" : [\"D1A2\", \"D3A3\", \"D4A7\"],\n",
    "        \"Tat6\" : [\"D1A1\", \"D2A3\"],\n",
    "        \"Tat4\" : [\"D1A1\", \"D2A2\", \"D3A3\", \"D4A7\"],\n",
    "        \"Tat9\" : [\"D1cA3\", \"D4A7\"],\n",
    "        \"Tat7\" : [\"D1A2\", \"D3A3\"],\n",
    "    \n",
    "        \"Tat8\" : [\"D1A1\", \"D2A2\", \"D3A3\"],\n",
    "        \"Tat10\" : [\"D1cA1b\", \"D2A3\", \"D4A7\"],\n",
    "        \"Tat11\" : [\"D1cA1b\", \"D2A2\", \"D3A3\", \"D4A7\"],\n",
    "        \"Tat12\" : [\"D1cA3\"],\n",
    "        \"Tat13\" : [\"D1A1b\", \"D2A3\", \"D4A7\"],\n",
    "        \"Tat15\" : [\"D1cA1\", \"D2A3\", \"D4A7\"],\n",
    "        \"Tat16\" : [\"D1A1c\", \"D2A3\", \"D4A7\"],\n",
    "        \"TTat1\" : [\"D1A3\", \"D4aA7\"],\n",
    "        \"TTat2\" : [\"D1A1\", \"D2A3\", \"D4aA7\"],\n",
    "\n",
    "        \"Rev2\" : [\"D1A4a\", \"D4A7\"],\n",
    "        \"Rev1\" : [\"D1A4b\", \"D4A7\"],\n",
    "        \"Rev7\" : [\"D1A2\", \"D3A4b\", \"D4A7\"],\n",
    "        \"Rev8\" : [\"D1A2\", \"D3A4a\", \"D4A7\"],\n",
    "        \"Rev13\" : [\"D1cA4b\", \"D4A7\"],\n",
    "        \"Rev3\" : [\"D1A4c\", \"D4A7\"],\n",
    "        \"Rev4\" : [\"D1A1\", \"D2A4b\", \"D4A7\"],\n",
    "        \"Rev5\" : [\"D1A1\", \"D2A4a\", \"D4A7\"],\n",
    "        \"Rev11\" : [\"D1A1\", \"D2A2\", \"D3A4a\", \"D4A7\"],\n",
    "        \"Rev10\" : [\"D1A1\", \"D2A2\", \"D3A4b\", \"D4A7\"],\n",
    "        \"Rev9\" : [\"D1A2\", \"D3A4c\", \"D4A7\"],\n",
    "        \n",
    "        \"Rev6\" : [\"D1A1\", \"D2A4c\", \"D4A7\"],\n",
    "        \"Rev12\" : [\"D1A1\", \"D2A2\", \"D3A4c\", \"D4A7\"],\n",
    "        \"Rev4b\" : [\"D1A4a\", \"D4aA7\"],\n",
    "        \"TRev\" : [\"D1A4a\", \"D4aA7\"],\n",
    "    \n",
    "        \"Nef2\" : [\"D1A5\", \"D4A7\"],\n",
    "        \"Nef4\" : [\"D1A2\", \"D3A5\", \"D4A7\"],\n",
    "        \"Nef3\" : [\"D1A1\", \"D2A5\", \"D4A7\"],\n",
    "        \"Nef5\" : [\"D1A1\", \"D2A2\", \"D3A5\", \"D4A7\"],\n",
    "        \"Nef1\" : [\"D1A7\"],\n",
    "        \"Nef13\" : [\"D1A5\", \"D4aA7\"],\n",
    "        \"Nef9\" : [\"D1A5a\", \"D4A7\"],\n",
    "        \"Nef14\" : [\"D1A2\", \"D3A5\", \"D4aA7\"],\n",
    "        \"Nef12\" : [\"D1cA5\", \"D4A7\"],\n",
    "        \"Nef16\" : [\"D1cA5c\", \"D4A7\"],\n",
    "        \"Nef17\" : [\"D1A5\", \"D4bA7\"],\n",
    "    \n",
    "        \"Nef11\" : [\"D1A2\", \"D3A5a\", \"D4A7\"],\n",
    "        \"Nef15\" : [\"D1cA1b\", \"D2A5\", \"D4A7\"],\n",
    "        \"Nef18\" : [\"D1A5\", \"D4A7c\"],\n",
    "        \"Nef19\" : [\"D1A1\", \"D2A5a\", \"D4A7\"],\n",
    "        \"Nef20\" : [\"D1A1b\", \"D2A5\", \"D4A7\"],\n",
    "        \"Nef21\" : [\"D1A1\", \"D2A5\", \"D4aA7\"],\n",
    "        \"Nef22\" : [\"D1cA1\", \"D2A5\", \"D4A7\"],\n",
    "        \"Nef23\" : [\"D1A1\", \"D2A5b\", \"D4A7\"],\n",
    "        \"Nef25\" : [\"D1A1\", \"D2dA7\"],\n",
    "\n",
    "        \"sORF1\" : [\"D1A1\", \"D2bA3\", \"D4A7\"],\n",
    "        \"sORF3\" : [\"D1A1\", \"D2bA2\", \"D3A3\", \"D4A7\"]\n",
    "\n",
    "    }#\"2LTR\" : [\"\"]\n",
    "\n",
    "#sorted(isoforms.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "7e9b50ec-6384-4725-ad66-dc09719aced2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'1': 290,\n",
       " '1c': 294,\n",
       " '1d': 340,\n",
       " '2': 4509,\n",
       " '2c': 4520,\n",
       " '2d': 4553,\n",
       " '2b': 4605,\n",
       " '3': 5010,\n",
       " '4': 5591,\n",
       " '4a': 5595,\n",
       " '4b': 5598}"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "donor_sites"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "b0cfc640-caf9-4e87-99ec-bdf936a0794c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'1': 4460,\n",
       " '1b': 4466,\n",
       " '1c': 4469,\n",
       " '2': 4937,\n",
       " '3': 5324,\n",
       " '4c': 5483,\n",
       " '4a': 5501,\n",
       " '4b': 5507,\n",
       " '5': 5523,\n",
       " '5a': 5527,\n",
       " '5b': 5530,\n",
       " '5c': 5539,\n",
       " '7c': 7892,\n",
       " '7': 7916}"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "acceptor_sites"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "44908a8f-c61c-4378-a4bf-e7b6689857ec",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "88"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(isoforms)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "638163eb-a748-4601-a8a6-805d1254b0ff",
   "metadata": {},
   "source": [
    "# Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c3b5bb84-8dae-4de8-b62c-491ac188fd12",
   "metadata": {},
   "outputs": [],
   "source": [
    "import re\n",
    "def get_splice_boundaries(string):\n",
    "    m = re.match(\"D(.*)A(.*)\", string)\n",
    "    return donor_sites[m[1]], acceptor_sites[m[2]]\n",
    "\n",
    "#fix start and end of exons, translate dictionaries to lists\n",
    "def get_exons(exon_list):\n",
    "\n",
    "    start=1\n",
    "    end=9000\n",
    "\n",
    "    exons = [[start,end]]\n",
    "    for i, splicesite in enumerate(exon_list):\n",
    "        donor_site, acceptor_site = get_splice_boundaries(splicesite)\n",
    "        exons[i][1] = donor_site\n",
    "        exons.append([acceptor_site, end])\n",
    "\n",
    "    return exons\n",
    "\n",
    "def calc_isoform_length(exon_boundaries):\n",
    "    length = 0\n",
    "    for exon in exon_boundaries:\n",
    "        length += exon[1] - exon[0]\n",
    "    return length\n",
    "    \n",
    "#remove exons that come after PCR primer\n",
    "def split_exons_to_PCR_end(exon_boundaries, PCR_end):\n",
    "    cleaned_exon_list = []\n",
    "    for exon in exon_boundaries:\n",
    "        if exon[0] > PCR_end:\n",
    "            continue\n",
    "        elif exon[1] > PCR_end:\n",
    "            exon[1] = PCR_end\n",
    "        cleaned_exon_list.append(exon)\n",
    "    return cleaned_exon_list\n",
    "\n",
    "# remove exons that come after the RT primer\n",
    "def actual_exons(exons, actual_exon_boundaries):\n",
    "    return exons[:len(actual_exon_boundaries)-1]\n",
    "\n",
    "\n",
    "def has_D4A7(splice_sites):\n",
    "    return \"D4A7\" in splice_sites\n",
    "\n",
    "#check whether primer binding site is within an exon\n",
    "def check_primer_binding_site(exon_boundaries, primer_binding_site):\n",
    "    primer_binds = False\n",
    "    for exon in exon_boundaries:\n",
    "        if (exon[0] < primer_binding_site) & (exon[1] > primer_binding_site):\n",
    "            primer_binds = True\n",
    "    return primer_binds\n",
    "\n",
    "def get_seq_from_exons(refseq, exon_list):\n",
    "    seq = []\n",
    "    for exon in exon_list:\n",
    "        if exon[0]>1:\n",
    "            seq.extend(refseq[exon[0]-2:exon[1]-1])\n",
    "        else:\n",
    "            seq.extend(refseq[exon[0]-1:exon[1]-1])\n",
    "    return \"\".join(seq)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d89fbe40-436e-46b8-988e-f5477b488ee8",
   "metadata": {},
   "source": [
    "# Generate GTF file for RT-PCR US, PS, FS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce8e9542-5a58-4bdb-973a-a4d437aeebd4",
   "metadata": {},
   "outputs": [],
   "source": [
    "#############################################################################\n",
    "#hardcoded RT-PCR US or RT-PCR PS/FS end\n",
    "def get_PCR_end(isoform):\n",
    "    if isoform == \"unspliced1\":\n",
    "        return 537\n",
    "    else:\n",
    "        return 5545\n",
    "\n",
    "#add prefix to isoform if reverse transcribed by RT primer\n",
    "def add_splice_info_to_name(name, RT1, RT2, RT3, PCR):\n",
    "    if not PCR:\n",
    "        return \"noPCR_\" + name\n",
    "    \n",
    "    captured = False\n",
    "    if RT1:\n",
    "        name = \"RT1_\" + name\n",
    "        captured = True\n",
    "    if RT2:\n",
    "        name = \"RT2_\" + name\n",
    "        captured = True\n",
    "    if RT3:\n",
    "        name = \"RT3_\" + name\n",
    "        captured = True\n",
    "    \n",
    "    if not captured:\n",
    "        return \"noRT_\" + name\n",
    "    \n",
    "    return name\n",
    "\n",
    "def expected(name, RT1, RT2, RT3, PCR):\n",
    "    if not PCR:\n",
    "        return False\n",
    "    if RT1:\n",
    "        return True\n",
    "    if RT2:\n",
    "        return True\n",
    "    if RT3:\n",
    "        return True\n",
    "    \n",
    "    return False\n",
    "\n",
    "\n",
    "file_suffix = \"PCR1_2_3\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4af1525a-2920-4b12-8a78-4ba634b34328",
   "metadata": {},
   "source": [
    "# Generate dataframe and add information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d4d9071-a4ad-4cba-bca7-94eb73496a9c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "isoforms_df = pd.DataFrame(zip(isoforms.keys(), isoforms.values()), columns = [\"isoform\", \"splice_sites\"])\n",
    "\n",
    "isoforms_df[\"exon_boundaries\"] = isoforms_df[\"splice_sites\"].apply(lambda x: get_exons(x))\n",
    "isoforms_df[\"length\"] = isoforms_df[\"exon_boundaries\"].apply(lambda x: calc_isoform_length(x))\n",
    "isoforms_df[\"D4A7\"] = isoforms_df[\"splice_sites\"].apply(lambda x : has_D4A7(x))\n",
    "\n",
    "isoforms_df[\"RT1_primer\"] = isoforms_df[\"exon_boundaries\"].apply(lambda x: check_primer_binding_site(x,536)) #od\n",
    "isoforms_df[\"RT2_primer\"] = isoforms_df[\"exon_boundaries\"].apply(lambda x: check_primer_binding_site(x,5782)) #od\n",
    "isoforms_df[\"RT3_primer\"] = isoforms_df[\"D4A7\"]\n",
    "\n",
    "#in virion, may amplify full length\n",
    "isoforms_df[\"RT4_primer\"] = isoforms_df[\"exon_boundaries\"].apply(lambda x: check_primer_binding_site(x,7133)) #od1205\n",
    "\n",
    "isoforms_df[\"PCR_primer\"] = isoforms_df[\"exon_boundaries\"].apply(lambda x: check_primer_binding_site(x,5545))\n",
    "\n",
    "isoforms_df[\"PCR_end\"] = isoforms_df[\"isoform\"].apply(lambda x: get_PCR_end(x))\n",
    "\n",
    "isoforms_df.loc[isoforms_df[\"isoform\"] == \"unspliced1\", \"RT2_primer\"] = False\n",
    "isoforms_df.loc[isoforms_df[\"isoform\"] == \"unspliced2\", \"RT1_primer\"] = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b8300c3-effc-46ae-9aec-f100ba3f63be",
   "metadata": {},
   "outputs": [],
   "source": [
    "isoforms_df[\"actual_exon_boundaries\"] = isoforms_df.apply(lambda x: split_exons_to_PCR_end(x.exon_boundaries, x.PCR_end), axis=1)\n",
    "isoforms_df[\"actual_exon_boundaries_str\"] = isoforms_df[\"actual_exon_boundaries\"].apply(lambda x:str(x))\n",
    "isoforms_df[\"actual_splice_sites\"] = isoforms_df.apply(lambda x: actual_exons(x.splice_sites, x.actual_exon_boundaries), axis=1)\n",
    "isoforms_df[\"actual_length\"] = isoforms_df[\"actual_exon_boundaries\"].apply(lambda x: calc_isoform_length(x))\n",
    "\n",
    "isoforms_df[\"isoform_name\"] = isoforms_df.apply(lambda x: add_splice_info_to_name(x.isoform, x.RT1_primer, x.RT2_primer, x.RT3_primer, x.PCR_primer), axis=1)\n",
    "isoforms_df[\"expected\"] = isoforms_df.apply(lambda x: expected(x.isoform, x.RT1_primer, x.RT2_primer, x.RT3_primer, x.PCR_primer), axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76ab23a2-530b-4dfd-9245-6a4f07b1ffa6",
   "metadata": {},
   "outputs": [],
   "source": [
    "isoforms_df.to_excel(f\"references/all_isoforms_{file_suffix}.xlsx\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58622558-3de3-4216-80c7-846985ac3889",
   "metadata": {},
   "outputs": [],
   "source": [
    "filtered_isoforms_df = isoforms_df[isoforms_df[\"expected\"]]\n",
    "filtered_isoforms_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e1d334c-e288-4e94-a715-a8af58a6ab85",
   "metadata": {},
   "outputs": [],
   "source": [
    "#check if 2 or more isoforms have the exact same splice sites (because they only differ in a splice site not covered by the RT/PCR)\n",
    "collision_fix_dict = {}\n",
    "for splice_site, value in filtered_isoforms_df.value_counts(\"actual_exon_boundaries_str\").items():\n",
    "    if value > 1:\n",
    "        print(splice_site)\n",
    "\n",
    "        colliding_isoforms = filtered_isoforms_df[filtered_isoforms_df[\"actual_exon_boundaries_str\"] == splice_site][\"isoform_name\"].values\n",
    "        print(colliding_isoforms)\n",
    "        multispliced = filtered_isoforms_df[filtered_isoforms_df[\"actual_exon_boundaries_str\"] == splice_site][\"D4A7\"].values\n",
    "                \n",
    "        collision_fix_dict[splice_site] = \"_\".join(sorted(colliding_isoforms))\n",
    "    else:\n",
    "        collision_fix_dict[splice_site] = filtered_isoforms_df[filtered_isoforms_df[\"actual_exon_boundaries_str\"] == splice_site][\"isoform_name\"].values[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26530e69-6ee2-4362-88cb-6a0699e49a20",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print table with expected lengths for different RT-PCRs\n",
    "filtered_isoforms_df[filtered_isoforms_df[\"RT2_primer\"]][[\"isoform\", \"actual_length\"]].sort_values(\"actual_length\", ascending=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e90a36c5-fda8-40dd-889e-f29f0111f5ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "gtf_df = filtered_isoforms_df[[\"isoform_name\", \"actual_splice_sites\", \"actual_exon_boundaries\", \"actual_exon_boundaries_str\", \"actual_length\"]].copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "afe05d99-249f-4bd3-9832-6154f1ffce48",
   "metadata": {},
   "outputs": [],
   "source": [
    "gtf_df.drop_duplicates(subset =\"actual_exon_boundaries_str\", inplace = True)\n",
    "gtf_df[\"isoform_name\"] = gtf_df[\"actual_exon_boundaries_str\"].apply(lambda x:collision_fix_dict[x])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7bfdd3ee-f978-490f-96ae-920275571d31",
   "metadata": {},
   "outputs": [],
   "source": [
    "gtf_df.to_pickle(f\"references/gtf_dataframe_{file_suffix}.pickle\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ce5b55b-6173-4939-8990-c56104898b5a",
   "metadata": {},
   "outputs": [],
   "source": [
    "gtf_df.to_excel(f\"references/actual_isoforms_{file_suffix}.xlsx\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7af5aea7-c41e-4ad4-bd40-4b899b3e34e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "gtf_outfile = f\"references/HIVNL43_all_iso_{file_suffix}.gtf\"\n",
    "\n",
    "with open(gtf_outfile, \"w\") as outfile:\n",
    "\n",
    "    for _, transcript in gtf_df.sort_values(\"actual_exon_boundaries_str\", ascending=False).iterrows():\n",
    "        transcript_start = transcript.actual_exon_boundaries[0][0]\n",
    "        transcript_end = transcript.actual_exon_boundaries[-1][1]\n",
    "\n",
    "        #gff format:\n",
    "        #seqname\\tsource\\tfeature\\tstart\\tend\\tscore\\tstrand\\tframe\\tattribute\n",
    "        #attributes:\n",
    "        #gene_id \"\"; gene_name \"\"; transcript_id \"\"; transcript_type \"\"; exon_number \"\";\n",
    "        \n",
    "        seqname, source, feature, score, strand, frame = \"p12\", \"PMID32807178\", \"transcript\", \"1000\", \"+\", \".\"\n",
    "        start, end = transcript_start, transcript_end\n",
    "        attributes = f'gene_id \\\"HIVNL43\\\"; gene_name \\\"HIVNL43\\\"; transcript_id \\\"{transcript.isoform_name}\\\"; transcript_type \\\"protein_coding\\\"'\n",
    "        \n",
    "        print(seqname, source, feature, start, end, score, strand, frame, attributes, sep=\"\\t\", end=\"\\n\", file=outfile)\n",
    "\n",
    "        for exon_number, exon in enumerate(transcript.actual_exon_boundaries):\n",
    "            exon_start = exon[0]\n",
    "            exon_end = exon[1]\n",
    "            start, end = exon_start, exon_end\n",
    "            feature = \"exon\"\n",
    "            attributes = f'gene_id \\\"HIVNL43\\\"; gene_name \\\"HIVNL43\\\"; transcript_id \\\"{transcript.isoform_name}\\\"; transcript_type \\\"protein_coding\\\"; exon_number \\\"{exon_number+1}\\\"'\n",
    "            print(seqname, source, feature, start, end, score, strand, frame, attributes, sep=\"\\t\", end=\"\\n\", file=outfile)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a788a035-e730-4065-b01e-18cfacb1071d",
   "metadata": {},
   "source": [
    "## Create fasta transcripts for spliced isoforms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6485b74-afb5-4322-a8fb-3b82fae4ff98",
   "metadata": {},
   "outputs": [],
   "source": [
    "gtf_df[\"isoform_name\"].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fbabda03-cc1d-4e4d-9396-fe74f312a1da",
   "metadata": {},
   "outputs": [],
   "source": [
    "from Bio import SeqIO\n",
    "reference_file = \"references/p12_HIVNL43.fa\"\n",
    "refseq = list(SeqIO.read(reference_file, \"fasta\").seq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb15f063-1bd2-4b09-a242-0fcdb879a248",
   "metadata": {},
   "outputs": [],
   "source": [
    "from Bio.SeqRecord import SeqRecord\n",
    "from Bio.Seq import Seq\n",
    "\n",
    "individual_isoforms_folder = f\"references/transcripts_{file_suffix}/\"\n",
    "os.makedirs(individual_isoforms_folder, exist_ok = True)\n",
    "\n",
    "for _, transcript in gtf_df.sort_values(\"actual_exon_boundaries_str\", ascending=False).iterrows():\n",
    "    sequence = get_seq_from_exons(refseq, transcript[\"actual_exon_boundaries\"])\n",
    "    name = transcript[\"isoform_name\"]\n",
    "    record = SeqRecord(\n",
    "    Seq(sequence),\n",
    "    id = name,\n",
    "    name = name\n",
    "    )\n",
    "    SeqIO.write(record, f'{individual_isoforms_folder}/{name}.fa', \"fasta\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4d6ac89-36d7-49b4-ab1d-8d01d9fac09a",
   "metadata": {},
   "source": [
    "## Create mask file for primers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "85164cd1-3f8d-4777-81df-08bb4899eccd",
   "metadata": {},
   "outputs": [],
   "source": [
    "pcr1_fw = \"GGTCTCTCTGGTTAGACCAGATCTGAG\" \n",
    "pcr1_re = \"GACAAATACTGGGACAGCTACAACCATC\"\n",
    "pcr23_re = \"GAAGCGGAGACAGCGACGAA\" \n",
    "pcr5_1_re = \"GAAGTAATTCCAGCAGAGACAGGG\"\n",
    "pcr5_2_fw = \"GAAGTAATTCCAGCAGAGACAGGG\"\n",
    "pcr5_2_re = \"CAATAAAGCTTGCCTTGAGTGCTTC\"\n",
    "\n",
    "with open(f\"references/mask_primers_{file_suffix}.csv\", \"w\") as outfile:\n",
    "    for _, row in gtf_df.iterrows():\n",
    "        if row.isoform_name == \"RT1_unspliced1\":\n",
    "            print(row.isoform_name, pcr1_fw, pcr1_re, sep=\";\", file=outfile)\n",
    "        else:\n",
    "            print(row.isoform_name, pcr1_fw, pcr23_re, sep=\";\", file=outfile)\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79c26bf3-0fc2-432b-9705-28d395ed3997",
   "metadata": {},
   "source": [
    "# Create gtf, primer mask and transcript sequences for RT-PCR FL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d486e58-78ee-4747-80cd-4a60c729dbcb",
   "metadata": {},
   "outputs": [],
   "source": [
    "#############################################################################\n",
    "#hardcoded PCR5.1 and 5.2 end\n",
    "def get_PCR_end(isoform):\n",
    "\n",
    "    return 7904 #pcr5_2 od1265\n",
    "\n",
    "#add prefix to isoform if reverse transcribed by RT primer\n",
    "def add_splice_info_to_name(name, PCR5_1_primer_re, PCR5_2_primer_fw, PCR5_2_primer_re):\n",
    "\n",
    "    captured = False\n",
    "    if PCR5_1_primer_re:\n",
    "        name = \"RT5_1_\" + name\n",
    "        captured = True\n",
    "    if PCR5_2_primer_fw & PCR5_2_primer_re:\n",
    "        name = \"RT5_2_\" + name\n",
    "        captured = True\n",
    "    if not captured:\n",
    "        return \"noRT_\" + name\n",
    "    \n",
    "    return name\n",
    "\n",
    "def expected(name, PCR5_1_primer_re, PCR5_2_primer_fw, PCR5_2_primer_re):\n",
    "\n",
    "    if PCR5_1_primer_re:\n",
    "        return True\n",
    "    if PCR5_2_primer_fw & PCR5_2_primer_re:\n",
    "        return True\n",
    "    return False\n",
    "\n",
    "file_suffix = \"PCR5\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35ef415c-ba09-4659-b211-c43703ee9626",
   "metadata": {},
   "source": [
    "# Generate dataframe and add information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1be0d915-9110-456f-b2d9-453d58027b7e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "isoforms_df = pd.DataFrame(zip(isoforms.keys(), isoforms.values()), columns = [\"isoform\", \"splice_sites\"])\n",
    "\n",
    "isoforms_df[\"exon_boundaries\"] = isoforms_df[\"splice_sites\"].apply(lambda x: get_exons(x))\n",
    "isoforms_df[\"length\"] = isoforms_df[\"exon_boundaries\"].apply(lambda x: calc_isoform_length(x))\n",
    "isoforms_df[\"D4A7\"] = isoforms_df[\"splice_sites\"].apply(lambda x : has_D4A7(x))\n",
    "\n",
    "#reverse transcribed with oligo-dT, so all \n",
    "\n",
    "isoforms_df[\"PCR5_1_primer_re\"] = isoforms_df[\"exon_boundaries\"].apply(lambda x: check_primer_binding_site(x,4058)) #od1267\n",
    "isoforms_df[\"PCR5_2_primer_fw\"] = isoforms_df[\"exon_boundaries\"].apply(lambda x: check_primer_binding_site(x,4034)) #od1268\n",
    "isoforms_df[\"PCR5_2_primer_re\"] = isoforms_df[\"exon_boundaries\"].apply(lambda x: check_primer_binding_site(x,7904)) #od1265\n",
    "\n",
    "isoforms_df[\"PCR_end\"] = isoforms_df[\"isoform\"].apply(lambda x: get_PCR_end(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "786f506b-6d69-435c-9b22-eb8fec8ce26d",
   "metadata": {},
   "outputs": [],
   "source": [
    "isoforms_df[\"actual_exon_boundaries\"] = isoforms_df.apply(lambda x: split_exons_to_PCR_end(x.exon_boundaries, x.PCR_end), axis=1)\n",
    "isoforms_df[\"actual_exon_boundaries_str\"] = isoforms_df[\"actual_exon_boundaries\"].apply(lambda x:str(x))\n",
    "isoforms_df[\"actual_splice_sites\"] = isoforms_df.apply(lambda x: actual_exons(x.splice_sites, x.actual_exon_boundaries), axis=1)\n",
    "isoforms_df[\"actual_length\"] = isoforms_df[\"actual_exon_boundaries\"].apply(lambda x: calc_isoform_length(x))\n",
    "\n",
    "isoforms_df.drop(isoforms_df[isoforms_df[\"isoform\"] == \"unspliced2\"].index, inplace=True)\n",
    "                 \n",
    "isoforms_df[\"isoform_name\"] = isoforms_df.apply(lambda x: add_splice_info_to_name(x.isoform, x.PCR5_1_primer_re, x.PCR5_2_primer_fw, x.PCR5_2_primer_re), axis=1)\n",
    "isoforms_df[\"expected\"] = isoforms_df.apply(lambda x: expected(x.isoform, x.PCR5_1_primer_re, x.PCR5_2_primer_fw, x.PCR5_2_primer_re), axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15e956a9-41e2-45f4-8680-48e3cb46c67e",
   "metadata": {},
   "outputs": [],
   "source": [
    "filtered_isoforms_df = isoforms_df[isoforms_df[\"expected\"]]\n",
    "filtered_isoforms_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d2485c9-38df-42e2-aad6-32fc71b4083d",
   "metadata": {},
   "outputs": [],
   "source": [
    "filtered_isoforms_df.value_counts(\"actual_exon_boundaries_str\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "985baa9f-c625-4cef-9876-5d39ed36919b",
   "metadata": {},
   "outputs": [],
   "source": [
    "#check if 2 or more isoforms have the exact same splice sites (because they only differ in a splice site not covered by the RT/PCR)\n",
    "collision_fix_dict = {}\n",
    "for splice_site, value in filtered_isoforms_df.value_counts(\"actual_exon_boundaries_str\").items():\n",
    "    if value > 1:\n",
    "        print(splice_site)\n",
    "\n",
    "        colliding_isoforms = filtered_isoforms_df[filtered_isoforms_df[\"actual_exon_boundaries_str\"] == splice_site][\"isoform_name\"].values\n",
    "        print(colliding_isoforms)\n",
    "        multispliced = filtered_isoforms_df[filtered_isoforms_df[\"actual_exon_boundaries_str\"] == splice_site][\"D4A7\"].values\n",
    "                \n",
    "        collision_fix_dict[splice_site] = \"_\".join(sorted(colliding_isoforms))\n",
    "    else:\n",
    "        collision_fix_dict[splice_site] = filtered_isoforms_df[filtered_isoforms_df[\"actual_exon_boundaries_str\"] == splice_site][\"isoform_name\"].values[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9aa3d506-70c5-4451-9988-8bf43f18767a",
   "metadata": {},
   "outputs": [],
   "source": [
    "gtf_df = filtered_isoforms_df[[\"isoform_name\", \"actual_splice_sites\", \"actual_exon_boundaries\", \"actual_exon_boundaries_str\", \"actual_length\"]].copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1148765d-8c43-498d-b5f2-97ff330d9bcd",
   "metadata": {},
   "outputs": [],
   "source": [
    "gtf_df.drop_duplicates(subset =\"actual_exon_boundaries_str\", inplace = True)\n",
    "gtf_df[\"isoform_name\"] = gtf_df[\"actual_exon_boundaries_str\"].apply(lambda x:collision_fix_dict[x])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "49f70904-6106-4b78-9d1b-367cde6c14f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "gtf_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e4dbf3b-08a0-45e7-a074-f8b904a57262",
   "metadata": {},
   "outputs": [],
   "source": [
    "gtf_df.to_pickle(f\"references/gtf_dataframe_{file_suffix}.pickle\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97369707-62fc-400a-8d4b-72b55e671e0e",
   "metadata": {},
   "outputs": [],
   "source": [
    "gtf_outfile = f\"references/HIVNL43_PCR5_{file_suffix}.gtf\"\n",
    "\n",
    "with open(gtf_outfile, \"w\") as outfile:\n",
    "\n",
    "    for _, transcript in gtf_df.sort_values(\"actual_exon_boundaries_str\", ascending=False).iterrows():\n",
    "        transcript_start = transcript.actual_exon_boundaries[0][0]\n",
    "        transcript_end = transcript.actual_exon_boundaries[-1][1]\n",
    "\n",
    "        #gff format:\n",
    "        #seqname\\tsource\\tfeature\\tstart\\tend\\tscore\\tstrand\\tframe\\tattribute\n",
    "        #attributes:\n",
    "        #gene_id \"\"; gene_name \"\"; transcript_id \"\"; transcript_type \"\"; exon_number \"\";\n",
    "        \n",
    "        seqname, source, feature, score, strand, frame = \"p12\", \"PMID32807178\", \"transcript\", \"1000\", \"+\", \".\"\n",
    "        start, end = transcript_start, transcript_end\n",
    "        attributes = f'gene_id \\\"HIVNL43\\\"; gene_name \\\"HIVNL43\\\"; transcript_id \\\"{transcript.isoform_name}\\\"; transcript_type \\\"protein_coding\\\"'\n",
    "        \n",
    "        print(seqname, source, feature, start, end, score, strand, frame, attributes, sep=\"\\t\", end=\"\\n\", file=outfile)\n",
    "\n",
    "        for exon_number, exon in enumerate(transcript.actual_exon_boundaries):\n",
    "            exon_start = exon[0]\n",
    "            exon_end = exon[1]\n",
    "            start, end = exon_start, exon_end\n",
    "            feature = \"exon\"\n",
    "            attributes = f'gene_id \\\"HIVNL43\\\"; gene_name \\\"HIVNL43\\\"; transcript_id \\\"{transcript.isoform_name}\\\"; transcript_type \\\"protein_coding\\\"; exon_number \\\"{exon_number+1}\\\"'\n",
    "            print(seqname, source, feature, start, end, score, strand, frame, attributes, sep=\"\\t\", end=\"\\n\", file=outfile)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1f75d9a-e24e-4f43-b83c-5e30c451021a",
   "metadata": {},
   "source": [
    "## Create fasta transcripts for spliced isoforms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c714c0c4-c141-4ad7-83b0-6cf55ac587e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "gtf_df[\"isoform_name\"].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7eb2038-9314-4eed-94d4-4d51675d1a63",
   "metadata": {},
   "outputs": [],
   "source": [
    "from Bio import SeqIO\n",
    "reference_file = \"references/p12_HIVNL43.fa\"\n",
    "refseq = list(SeqIO.read(reference_file, \"fasta\").seq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3dc05b3a-ecfa-464d-8c90-f948e0cfa98f",
   "metadata": {},
   "outputs": [],
   "source": [
    "from Bio.SeqRecord import SeqRecord\n",
    "from Bio.Seq import Seq\n",
    "\n",
    "individual_isoforms_folder = f\"references/transcripts_{file_suffix}/\"\n",
    "os.makedirs(individual_isoforms_folder, exist_ok = True)\n",
    "\n",
    "for _, transcript in gtf_df.sort_values(\"actual_exon_boundaries_str\", ascending=False).iterrows():\n",
    "    sequence = get_seq_from_exons(refseq, transcript[\"actual_exon_boundaries\"])\n",
    "    name = transcript[\"isoform_name\"]\n",
    "    record = SeqRecord(\n",
    "    Seq(sequence),\n",
    "    id = name,\n",
    "    name = name\n",
    "    )\n",
    "    SeqIO.write(record, f'{individual_isoforms_folder}/{name}.fa', \"fasta\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cfe7bfec-db73-4617-b5f2-b9449ad8cd8d",
   "metadata": {},
   "source": [
    "## Create mask file for primers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "88e692c7-7e21-48b0-bd1d-f608ca958411",
   "metadata": {},
   "outputs": [],
   "source": [
    "pcr1_fw = \"GGTCTCTCTGGTTAGACCAGATCTGAG\"\n",
    "pcr1_re = \"GACAAATACTGGGACAGCTACAACCATC\" # in fw orientation, so revcomp of re primer\n",
    "pcr23_re = \"GAAGCGGAGACAGCGACGAA\"  # in fw orientation, so revcomp of re primer\n",
    "pcr5_1_re_2_fw = \"GAAGTAATTCCAGCAGAGACAGGG\"\n",
    "pcr5_2_re = \"CAATAAAGCTTGCCTTGAGTGCTTC\"  # in fw orientation, so revcomp of re primer\n",
    "\n",
    "with open(f\"references/mask_primers_{file_suffix}.csv\", \"w\") as outfile:\n",
    "    for _, row in gtf_df.iterrows():\n",
    "        print(row.isoform_name, pcr1_fw, pcr5_1_re_2_fw, pcr5_2_re, sep=\";\", file=outfile)\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee140bf0-4284-4d4a-86d5-e00b29beeda4",
   "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}