1.2_orf_selection_&_sitemapping.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Important\n",
    "\n",
    "`run_mafft.sh` in `SNPs_comparison` to make alignment before this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import json\n",
    "\n",
    "import pandas as pd\n",
    "from Bio import SeqIO, AlignIO\n",
    "from Bio.Align import MultipleSeqAlignment\n",
    "\n",
    "dataDir = \"Data\"\n",
    "snpCompareDir = \"SNPs_comparison\"\n",
    "snpDir = \"SNPs\"\n",
    "genomeTreeEXDir = \"Genome_tree_EX\"\n",
    "\n",
    "reference = \"EPI_ISL_402125\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "aligned = AlignIO.read(os.path.join(snpCompareDir, \"aligned.fasta\"), \"fasta\")\n",
    "annoRef = SeqIO.read(os.path.join(dataDir, \"reference.gb\"), \"gb\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "ref2aligned = {}\n",
    "\n",
    "for record in aligned:\n",
    "    if record.id == reference:\n",
    "        refIndex = -1\n",
    "        for alignedIndex, n in enumerate(record):\n",
    "            if n != '-':\n",
    "                refIndex += 1\n",
    "            ref2aligned[refIndex] = alignedIndex\n",
    "        break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "cdsAlignedIndex = []\n",
    "product_id = 0\n",
    "\n",
    "for f in annoRef.features:\n",
    "    if f.type == \"CDS\":\n",
    "        product_id += 1\n",
    "        for i in f.location.parts:\n",
    "            cdsAlignedIndex.append({\n",
    "                \"start\": ref2aligned[i.start],\n",
    "                \"end\": ref2aligned[i.end],\n",
    "                \"product_id\": product_id,\n",
    "                \"product_name\": \", \".join(f.qualifiers[\"product\"]),\n",
    "                \"ref_start\": i.start,\n",
    "                \"ref_end\": i.end\n",
    "            })"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alignedCDS = None\n",
    "cdsPos = 0\n",
    "cdsAnno = {}\n",
    "\n",
    "for loc in cdsAlignedIndex:\n",
    "    cdsSeq = aligned[:, loc[\"start\"]:loc[\"end\"]]\n",
    "    if alignedCDS is None:\n",
    "        alignedCDS = cdsSeq\n",
    "    else:\n",
    "        alignedCDS += cdsSeq\n",
    "\n",
    "AlignIO.write(alignedCDS, os.path.join(snpCompareDir, \"aligned_cds.fasta\"), \"fasta\")\n",
    "AlignIO.write(alignedCDS, os.path.join(genomeTreeEXDir, \"aligned_cds.fasta\"), \"fasta\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "ref2alignedCDS = {}\n",
    "ref2alignedPos = {}\n",
    "\n",
    "for record in alignedCDS:\n",
    "    if record.id == reference:\n",
    "        refIndex = -1\n",
    "        for alignedIndex, n in enumerate(record):\n",
    "            if n != '-':\n",
    "                refIndex += 1\n",
    "            ref2alignedCDS[refIndex] = alignedIndex\n",
    "            ref2alignedPos[refIndex + 1] = alignedIndex + 1\n",
    "        break\n",
    "\n",
    "with open(os.path.join(snpCompareDir, \"ref2aligned_cds.json\"), 'w') as f:\n",
    "    json.dump(ref2alignedPos, f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "136"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "snpPos = 8517\n",
    "pos = ref2alignedCDS[snpPos]\n",
    "seq = alignedCDS[:, pos - 15:pos + 12]\n",
    "\n",
    "outDir = os.path.join(snpCompareDir, str(snpPos))\n",
    "if not os.path.exists(outDir):\n",
    "    os.mkdir(outDir)\n",
    "    \n",
    "SeqIO.write(seq, os.path.join(outDir, \"nt{}.fasta\".format(snpPos)), \"fasta\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "nt = pd.DataFrame(seq)\n",
    "nt.columns = nt.columns.values + (snpPos - 14)\n",
    "nt.index = [record.id for record in seq]\n",
    "nt.to_csv(os.path.join(outDir, \"nt{}.csv\".format(snpPos)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "aaSeq = MultipleSeqAlignment([seq[0].translate(id=seq[0].id)])\n",
    "for s in seq[1:]:\n",
    "    aaSeq.add_sequence(s.id, str(s.translate().seq))\n",
    "\n",
    "SeqIO.write(aaSeq, os.path.join(outDir, \"aa{}.fasta\".format(snpPos)), \"fasta\")\n",
    "\n",
    "aa = pd.DataFrame(aaSeq)\n",
    "aa.columns = aa.columns.values + int((snpPos - 12) / 3)\n",
    "aa.index = [record.id for record in aaSeq]\n",
    "aa.to_csv(os.path.join(outDir, \"aa{}.csv\".format(snpPos)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "136"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "snpPos = 27641\n",
    "pos = ref2alignedCDS[snpPos]\n",
    "seq = alignedCDS[:, pos - 14:pos + 13]\n",
    "\n",
    "outDir = os.path.join(snpCompareDir, str(snpPos))\n",
    "if not os.path.exists(outDir):\n",
    "    os.mkdir(outDir)\n",
    "\n",
    "SeqIO.write(seq, os.path.join(outDir, \"nt{}.fasta\".format(snpPos)), \"fasta\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "nt = pd.DataFrame(seq)\n",
    "nt.columns = nt.columns.values + (snpPos - 13)\n",
    "nt.index = [record.id for record in seq]\n",
    "nt.to_csv(os.path.join(outDir, \"nt{}.csv\".format(snpPos)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "aaSeq = MultipleSeqAlignment([seq[0].translate(id=seq[0].id)])\n",
    "for s in seq[1:]:\n",
    "    aaSeq.add_sequence(s.id, str(s.translate().seq))\n",
    "\n",
    "SeqIO.write(aaSeq, os.path.join(outDir, \"aa{}.fasta\".format(snpPos)), \"fasta\")\n",
    "\n",
    "aa = pd.DataFrame(aaSeq)\n",
    "aa.columns = aa.columns.values + int((snpPos - 11) / 3)\n",
    "aa.index = [record.id for record in aaSeq]\n",
    "aa.to_csv(os.path.join(outDir, \"aa{}.csv\".format(snpPos)))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}