DiscobaMMSeqs2.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "DiscobaMMSeqs2.ipynb",
      "private_outputs": true,
      "provenance": [],
      "collapsed_sections": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/zephyris/discoba_alphafold/blob/main/DiscobaMMSeqs2.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ZOa_RGpZoD_0",
        "cellView": "form"
      },
      "source": [
        "#@title Generate a3m alignment using MMSeqs2 and a diverse set of Discoba species\n",
        "#@markdown Ideal for generating alignment for AlphaFold2 for kinetoplastid species like _Leishmania_ and _Trypanosoma_\n",
        "\n",
        "#@markdown Enter the query name and sequence then press <kbd>Ctrl</kbd> + <kbd>F9</kbd> or select `Runtime` > `Run All`\n",
        "from google.colab import files\n",
        "import re\n",
        "import hashlib\n",
        "\n",
        "query_name = 'Tb927.8.4970_PFR2' #@param {type:\"string\"}\n",
        "query_sequence = 'MLGTVDAIDYDGDRLHKVVLRFPAVRSGESEIVKEVWPCERIGQGSFGTVYRAVSSDYPRLALKISTGKSTRLRQELDVLSRVCTKGRLLLPRFEFGALNKTADLIVIGMELCVPSTLHDLLLSTRITSEAEMLFMAHQAVQAVSYVHAEGCIHRDIKLQNFVFDLDGNLKLIDFGLACNSLKPPAGDVVAGTVSFMSPEMAHNALHKDRRVSVGVAADVWSLGIVLFSIFTQRNPYPAPETPAPAAGSTPGGAGAAGVTGRGDITHGAEGEKGNDLSQQHRMNERLLRRVAAGDWQWPVGVTVSQDLKQLVNSILVVNPEERPSVSTILENKLWNLRRRYPPAAVAAFLGVQDDFLLSHDEAHLMRAVEERSAGVAASLRNSRLHSPASASNEDNGETDAQHSSSGSARNGGLNTSGATTSPARSSLKVVQRCGEGGIDGAVTVQVYDVRASTRKRSKPIREISVVMAEETAKTRRSKSARRATGAVSAPSSRVVSRAASTEYSRRIAPPAGARLQSSAAHTCANSGDDEGEAEAEAVNRGTSTSQRHSRGMSPVRLQDALETAGSVAAGQRDSVGHPKALVEPSCATPPLLTSGKQQPLSEHMRPDIQRSGSVELLEDAEAPTEASAMPTSHKRAASSGKKRRDASLRQPSSLILKGSTRDLSADAPRSTTTAASTQASTSLLASRTTLPLSAVNPSPSSSRQASLRRQASASAAAVSSAQGCAGHRGSSPVMKRAQRVALELGLDVIWHDEADHRRALSAMLLIEHAWLLASFRLTIEEDQERYSITWLAEEQEKSAAHPHRFKEVMQVMSKKYQYGFVCDMCDYEFLPTGPGEKDLHFFHCPCGRDLCPDCYTAYQRQCTCSCCRAVHSNSCVLREHLLLTGGTQYYSGSRKTNAAARADAVRGSFQAAASLNEEAESGDEASAPPEPPRRRGRPPKQDKNRSAVKQKGSRAAKDSSRRRRGAQDTLDVSVDDAHEVEQINLPRISIAAMQQQEERSSNGSHRGGGTAAVGVAPRPQRPEDVEVKQRPVESVPEGPWRPFARFKKDRRDEVAQQPTPEERDALLNGEWIRHFYLFPQAEPERVAASGTWAEGEEEPYAFVYHAQPGRTGAIFLTSDFPMHSAVFSMLERQFFVVNQVDTVEGVDSTRATSLLKAKGHPELRIAFHALQDIVAYDTNMMKQQRTPGTVSVYQAPRSAYSCNGEPFLYVRWFRFNENRTLSAFLLSNGAVQVFVNNEYELRWFDESRKFLIRYNGVCELVDDGTFALAPGINHLLYDSFDA' #@param {type:\"string\"}\n",
        "\n",
        "def add_hash(x,y):\n",
        "  return x+\"_\"+hashlib.sha1(y.encode()).hexdigest()[:5]\n",
        "\n",
        "query_name = \"\".join(query_name.split())\n",
        "query_name = re.sub(r'\\W+', '', query_name)\n",
        "query_name = add_hash(query_name, query_sequence)\n",
        "\n",
        "with open(f\"{query_name}.fasta\", \"w\") as text_file:\n",
        "    text_file.write(\">1\\n%s\" % query_sequence)\n",
        "\n",
        "#@markdown Searches take several minutes. The first run will set up the MMSeqs2 software and download a large sequence database so will take a little longer.\n",
        "\n",
        "#@markdown A full description of this resource is available in [this paper](https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0259871).\n",
        "#@markdown If you make use of this resource please cite: Wheeler RJ. \"A resource for improved predictions of Trypanosoma and Leishmania protein three-dimensional structure\" _PLoS One_, doi:&nbsp;[10.1371/journal.pone.0259871](https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0259871) (2021).\n",
        "\n",
        "#@markdown You should also cite MMSeqs2, information about citing MMSeqs2 is available on their [Github page](https://github.com/soedinglab/MMseqs2)."
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ixyncLlwkK6B",
        "cellView": "form"
      },
      "source": [
        "#@title Install MMseqs2, fetch and build Discoba database\n",
        "#@markdown Downloads a database database which comprehensively gathers Discoba protein sequence databases.\n",
        "\n",
        "#@markdown Many of these sequences are not in the UniRef100 database, which is commonly used by online database search/alignment tools.\n",
        "\n",
        "%%bash -s\n",
        "#Install MMSEQS2, for search and alignment\n",
        "if [ ! -f MMSEQSLOC_READY ]; then\n",
        "  curl https://mmseqs.com/latest/mmseqs-linux-avx2.tar.gz -s -L -o mmseqs-linux-avx2.tar.gz\n",
        "  tar -xf mmseqs-linux-avx2.tar.gz\n",
        "  touch MMSEQSLOC_READY\n",
        "fi\n",
        "\n",
        "#Install CD-HIT, for reducing complexity of the database\n",
        "#if [ ! -f CDHIT_READY ]; then\n",
        "#  git clone https://github.com/weizhongli/cdhit\n",
        "#  cd cdhit\n",
        "#  make\n",
        "#  cd ..\n",
        "#  touch CDHIT_READY\n",
        "#fi\n",
        "\n",
        "#Download the custom Discoba database\n",
        "if [ ! -f DISCOBA_READY ]; then\n",
        "  if [ -d discoba ]; then\n",
        "    rm -r discoba\n",
        "  fi\n",
        "  mkdir discoba\n",
        "  cd discoba\n",
        "    curl https://zenodo.org/record/5682928/files/discobaStats.txt?download=1\n",
        "    curl https://zenodo.org/record/5682928/files/discoba.fasta.gz?download=1 -s -L -o discoba.fasta.gz\n",
        "    gzip -d discoba.fasta.gz\n",
        "    #cdhit/cd-hit-est -i discoba.full.fasta -o discoba.fasta\n",
        "    ../mmseqs/bin/mmseqs createdb discoba.fasta discoba -v 0\n",
        "  cd ..\n",
        "  touch DISCOBA_READY\n",
        "fi"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "syldUPAooGrD",
        "cellView": "form"
      },
      "source": [
        "#@title Run MMseqs2 to get Discoba MSA\n",
        "#@markdown It is normal for this to take several minutes.\n",
        "%%bash -s $query_name\n",
        "MMSEQS=../mmseqs/bin/mmseqs\n",
        "DISCOBA=../discoba/discoba\n",
        "NAME=$1\n",
        "if [ ! -f $NAME.mmseqs ]; then\n",
        "  if [ -d mmseqs_$NAME ]; then\n",
        "    rm -r mmseqs_$NAME\n",
        "  fi\n",
        "  mkdir mmseqs_$NAME\n",
        "  cd mmseqs_$NAME\n",
        "  cp ../$NAME.fasta query.fasta\n",
        "  $MMSEQS createdb query.fasta query -v 1\n",
        "  echo \"Searching database\"\n",
        "  $MMSEQS search query $DISCOBA search tmp --num-iterations 2 -v 1\n",
        "  echo \"Aligning hits\"\n",
        "  $MMSEQS align query $DISCOBA search align -a -v 1\n",
        "  echo \"Formatting result\"\n",
        "  $MMSEQS convertalis query $DISCOBA align query.tab --format-output target,qlen,qstart,qend,tstart,tend,tseq,cigar,taln -v 0\n",
        "  touch $NAME.mmseqs\n",
        "fi"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Wk9qbiHzpJfE",
        "cellView": "form"
      },
      "source": [
        "#@title Reformat MMseq2 alignment table to .a3m\n",
        "mmseq_tab=open(\"mmseqs_\"+query_name+\"/query.tab\", \"r\")\n",
        "mmseq_data=mmseq_tab.readlines()\n",
        "mmseq_out=open(query_name+\".a3m\", \"w\")\n",
        "mmseq_count=0\n",
        "for mmseq_line in mmseq_data:\n",
        "  mmseq_line=mmseq_line.replace(\"\\r\", \"\").replace(\"\\n\", \"\").split(\"\\t\")\n",
        "  mmseq_count+=1\n",
        "  #Parse cigar\n",
        "  mmseq_cigar=re.findall('(\\d+)([MDI])?', mmseq_line[7])\n",
        "  #Pad start of alignment\n",
        "  alignment_seq=\"-\"*(int(mmseq_line[2])-1)\n",
        "  alignment_index=int(mmseq_line[4])-1;\n",
        "  #Loop through cigar re-writing sequence\n",
        "  for cigar_entry in mmseq_cigar:\n",
        "    for sequence_index in range(0, int(cigar_entry[0])):\n",
        "      if cigar_entry[1]==\"M\":\n",
        "        alignment_seq+=mmseq_line[6][alignment_index:alignment_index+1].upper()\n",
        "        alignment_index+=1\n",
        "      elif cigar_entry[1]==\"D\":\n",
        "        alignment_seq+=mmseq_line[6][alignment_index:alignment_index+1].lower()\n",
        "        alignment_index+=1\n",
        "      elif cigar_entry[1]==\"I\":\n",
        "        alignment_seq+=\"-\"\n",
        "  #Pad end of alignment\n",
        "  alignment_seq=alignment_seq+(\"-\"*(int(mmseq_line[1])-int(mmseq_line[3])))\n",
        "  #Print result\n",
        "  mmseq_out.write(\">%s\\n\" % mmseq_line[0])\n",
        "  mmseq_out.write(\"%s\\n\" % alignment_seq)\n",
        "mmseq_out.close()\n",
        "print(\"%d sequences found\" % mmseq_count)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Ka0EGSc0F4rW",
        "cellView": "form"
      },
      "source": [
        "#@title Download the a3m file\n",
        "#@markdown This a3m file can be used as a standalone alignment for your favourite AlphaFold2 implementation.\n",
        "#@markdown It uses MMSeqs2 for alignment, as used in [this colab notebook](https://colab.research.google.com/github/sokrypton/ColabFold/blob/main/AlphaFold2.ipynb).\n",
        "\n",
        "#@markdown The content of this a3m can be appended to a a3m from searching a different database, so long as the same input sequence was used.\n",
        "files.download(f\"{query_name}.a3m\")"
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}