Bos_Project_Final.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bos taurus mitokondriumok kódoló mutációi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import re\n",
    "from math import floor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "infile=open('bos_sequence.fasta','r' )\n",
    "adat=infile.read()\n",
    "sorszerint=adat.split(\"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Teljes adatokat az egyes mitokondrium szekvenciák szerint tördeljük. (Címsorok már nem lesznek benne a tömbbben.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "299\n"
     ]
    }
   ],
   "source": [
    "AllMitoSeq=[]\n",
    "seq=\"\"\n",
    "elsosore=True\n",
    "for sor in sorszerint:\n",
    "    if len(sor)>0:\n",
    "        if sor[0]=='>':\n",
    "            if elsosore!=True:\n",
    "                AllMitoSeq.append(seq)\n",
    "            else:\n",
    "                elsosore=False\n",
    "            seq=\"\"\n",
    "        else:\n",
    "            seq+=sor\n",
    "AllMitoSeq.append(seq)\n",
    "    \n",
    "print(len(AllMitoSeq))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Leghosszabb és legrövidebb szekvenciák hosszának meghatározása."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "16337 16448\n"
     ]
    }
   ],
   "source": [
    "LenSeq=[]\n",
    "for seq in AllMitoSeq:\n",
    "    LenSeq.append(len(seq))\n",
    "    \n",
    "print(min(LenSeq),max(LenSeq))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A minimum és maximum közé eső (és megegyező) hosszúságú szekvenciákból hány darab van."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[16337, 10], [16338, 27], [16339, 87], [16340, 115], [16341, 12], [16342, 43], [16344, 1], [16356, 1], [16383, 1], [16426, 1], [16448, 1]]\n"
     ]
    }
   ],
   "source": [
    "LenStat=[]\n",
    "\n",
    "for i in np.arange(min(LenSeq),max(LenSeq)+1):\n",
    "    if LenSeq.count(i)>0:\n",
    "        LenStat.append([i,LenSeq.count(i)])\n",
    "\n",
    "print(LenStat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A leggyakoribb hosszúságú szekvenciák kiválogatása."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[16340, 115]\n"
     ]
    }
   ],
   "source": [
    "maxcount=0\n",
    "\n",
    "for i in range(len(LenStat)):\n",
    "    if LenStat[i][1]>maxcount:\n",
    "        maxcount=LenStat[i][1]\n",
    "        maxindex=i\n",
    "    \n",
    "print(LenStat[maxindex])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mégis a 16338 hosszúakat használtam, mert a [referencia genom](http://www.ncbi.nlm.nih.gov/nuccore/NC_006853.1) is ilyen hoszzú."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Mito=[]\n",
    "\n",
    "for sor in AllMitoSeq:\n",
    "    if len(sor)==(LenStat[1][0]):\n",
    "        Mito.append(sor)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Különbségi mátrix: Minden szekvenciát minden másik szekvenciával összehasonlítom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def compare(x1,x2):\n",
    "    k=0\n",
    "\n",
    "    for i in range(len(x1)):\n",
    "        if x1[i]!=x2[i]:\n",
    "            k=k+1\n",
    "    return k"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "matrixom=np.zeros((LenStat[maxindex][1],LenStat[maxindex][1]))\n",
    "\n",
    "for i in range(len(Mito)):\n",
    "    for j in range(len(Mito)):\n",
    "        matrixom[i,j]=compare(Mito[i],Mito[j])\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A ftp://ftp.ncbi.nlm.nih.gov/genomes/Bos_taurus/CHR_MT/bt_ref_Bos_taurus_UMD_3.1.1_chrMT.gbk.gz fileból kiszedem a géneket kódoló régiók kezdeti és végső nukleotidjának pozícióit ('gene'-el kezdődő sorok)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "reffile=open('bt_ref_Bos_taurus_UMD_3.1.1_chrMT.gbk', 'r')\n",
    "adat2=reffile.read()\n",
    "byline=adat2.split(\"\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "GeneSeq=[]\n",
    "\n",
    "for i in range(len(byline)):\n",
    "    if byline[i][0:9]=='     gene':\n",
    "        kulon=byline[i].split()\n",
    "        if kulon[1][0]!='c':\n",
    "            #A complementtel kezdődőt kidobtam\n",
    "            GeneSeq.append(kulon[1])\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[3101, 4266, 5687, 7374, 8129, 8290, 8970, 9823, 10239, 10529, 12109, 14514] [4056, 5307, 7231, 8057, 8329, 8970, 9750, 10168, 10535, 11906, 13929, 15653]\n"
     ]
    }
   ],
   "source": [
    "Start=[]\n",
    "End=[]\n",
    "\n",
    "for pos in GeneSeq:\n",
    "    startend=pos.split('..')\n",
    "    Start.append(int(startend[0]))\n",
    "    End.append(int(startend[1]))\n",
    "    \n",
    "print(Start,End)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mito-ban az adott, referencia genommal egyező hosszúságú mktokondrium szekvenciák vannak (26)  \n",
    "MitoGene-ben ezek vannak gének szerint szétválogatva  \n",
    "- Első index:egyedek\n",
    "- Második index:gének"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "MitoGene=[]\n",
    "for j in range(len(Mito)):\n",
    "    Individ=[]\n",
    "    for i in range(len(Start)):\n",
    "        Individ.append(Mito[j][Start[i]-1:End[i]-1])\n",
    "        if i == len(Start)-1:\n",
    "            MitoGene.append(Individ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "A referencia genomból is kiválogatom a géneket."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "infile=open('bos_ref_sequence.fasta','r' )\n",
    "adat=infile.read()\n",
    "sorszerint=adat.split(\"\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "RefSeq=\"\"\n",
    "\n",
    "for sor in sorszerint:\n",
    "    if len(sor)>0:\n",
    "        if sor[0]!='>':\n",
    "            RefSeq+=sor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "RefMitoGene=[]\n",
    "for i in range(len(Start)):\n",
    "    RefMitoGene.append(RefSeq[Start[i]-1:End[i]-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "És ezeket összehasonlítom az összes mitokondrium génszekvenciával."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "szekvmatrix=np.zeros((len(RefMitoGene),len(MitoGene)))\n",
    "\n",
    "for i in range(len(RefMitoGene)):\n",
    "    for j in range(len(MitoGene)):\n",
    "        szekvmatrix[i,j]=compare(RefMitoGene[i],MitoGene[j][i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Elejétől indulva aminosavakká átkonvertálom a szekvenciákat. és ezeket is összehasonlítom.  \n",
    "Nem minden gén osztható 3-al, de az egyszerűség kedvéért az elejétől indulva addig veszem, amíg van 3 bázis egymás után"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "map = {\"TTT\":\"Phe\", \"TTC\":\"Phe\", \"TTA\":\"Leu\", \"TTG\":\"Leu\",\n",
    "    \"TCT\":\"Ser\", \"TCC\":\"Ser\", \"TCA\":\"Ser\", \"TCG\":\"Ser\",\n",
    "    \"TAT\":\"Tyr\", \"TAC\":\"Tyr\", \"TAA\":\"STOP\", \"TAG\":\"STOP\",\n",
    "    \"TGT\":\"Cys\", \"TGC\":\"Cys\", \"TGA\":\"STOP\", \"TGG\":\"Trp\",\n",
    "    \"CTT\":\"Leu\", \"CTC\":\"Leu\", \"CTA\":\"Leu\", \"CTG\":\"Leu\",\n",
    "    \"CCT\":\"Pro\", \"CCC\":\"Pro\", \"CCA\":\"Pro\", \"CCG\":\"Pro\",\n",
    "    \"CAT\":\"His\", \"CAC\":\"His\", \"CAA\":\"Gln\", \"CAG\":\"Gln\",\n",
    "    \"CGT\":\"Arg\", \"CGC\":\"Arg\", \"CGA\":\"Arg\", \"CGG\":\"Arg\",\n",
    "    \"ATT\":\"Ile\", \"ATC\":\"Ile\", \"ATA\":\"Ile\", \"ATG\":\"Met\",\n",
    "    \"ACT\":\"Thr\", \"ACC\":\"Thr\", \"ACA\":\"Thr\", \"ACG\":\"Thr\",\n",
    "    \"AAT\":\"Asn\", \"AAC\":\"Asn\", \"AAA\":\"Lys\", \"AAG\":\"Lys\",\n",
    "    \"AGT\":\"Ser\", \"AGC\":\"Ser\", \"AGA\":\"Arg\", \"AGG\":\"Arg\",\n",
    "    \"GTT\":\"Val\", \"GTC\":\"Val\", \"GTA\":\"Val\", \"GTG\":\"Val\",\n",
    "    \"GCT\":\"Ala\", \"GCC\":\"Ala\", \"GCA\":\"Ala\", \"GCG\":\"Ala\",\n",
    "    \"GAT\":\"Asp\", \"GAC\":\"Asp\", \"GAA\":\"Glu\", \"GAG\":\"Glu\",\n",
    "    \"GGT\":\"Gly\", \"GGC\":\"Gly\", \"GGA\":\"Gly\", \"GGG\":\"Gly\",}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "RefAmino=[]\n",
    "\n",
    "for i in range(len(RefMitoGene)):\n",
    "    ByGene=[]\n",
    "    for j in range(0,int(floor((len(RefMitoGene[i])/3))*3),3):\n",
    "        kodon=\"\".join(RefMitoGene[i][j:j+3])\n",
    "        ByGene.append(map[kodon])\n",
    "        if j==(int(floor((len(RefMitoGene[i])/3))*3))-3:\n",
    "            RefAmino.append(ByGene)\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "AllAmino=[]\n",
    "\n",
    "for i in range(len(MitoGene)):\n",
    "    ByIndivid=[]\n",
    "    for x in range(len(MitoGene[i])):\n",
    "        ByGene=[]\n",
    "        for j in range(0,int((floor(len(MitoGene[i][x])/3)))*3,3):\n",
    "            kodon=\"\".join(MitoGene[i][x][j:j+3])\n",
    "            ByGene.append(map[kodon])\n",
    "            if j==(int(floor((len(MitoGene[i][x])/3))*3))-3:\n",
    "                ByIndivid.append(ByGene)\n",
    "                if x==len(MitoGene[i])-1:\n",
    "                    AllAmino.append(ByIndivid)\n",
    "            \n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Összehasonlítom a referencia aminosav génszekvenciákkal a minta aminosav génszekvanciákkal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "aminomatrix=np.zeros((len(RefAmino),len(AllAmino)))\n",
    "\n",
    "\n",
    "for i in range(len(RefAmino)):\n",
    "    for j in range(len(AllAmino)):\n",
    "        aminomatrix[i,j]=compare(RefAmino[i],AllAmino[j][i])\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A referencia szekvencia génjeit összehasonlítom a mintagénekkel, és ha a tripletek által kódolt aminosavak különböznek, kiírom a tripletet és annak helyét.  \n",
    "Az ujb-ben vannak az aminosav szekvenciában is jelen levő mutációk helyei."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "AllDiff = []\n",
    "for i in range(len(MitoGene)):\n",
    "    for g in range(len(RefMitoGene)):\n",
    "        lista=[]\n",
    "        for b in range(0,int(floor((len(RefMitoGene[g])/3))*3),3):\n",
    "            RefCodon=\"\".join(RefMitoGene[g][b:b+3])\n",
    "            AllCodon=\"\".join(MitoGene[i][g][b:b+3])\n",
    "            RefAA = map[RefCodon]\n",
    "            AllAA = map[AllCodon]\n",
    "            if RefAA!=AllAA:\n",
    "                for nn in range(3):\n",
    "                    if RefCodon[nn]!=AllCodon[nn]:\n",
    "                        ujb=b+nn\n",
    "                        AllDiff.append([i,g,b,ujb,RefCodon,AllCodon,RefAA,AllAA])\n",
    "                "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "infile=open('muthelyek.csv',\"w\")\n",
    "infile.write(\"Egyed,Gen,Tripletkezdet,Mutans bazishely,Referencia triplet,Minta triplet,Referencia aminosav,Minta aminosav \\n\")\n",
    "for i in AllDiff:\n",
    "    infile.write(str(i).replace('[','').replace(']','')+'\\n')\n",
    "infile.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}