BLAST-SprotVariants.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6bd639a6",
   "metadata": {},
   "source": [
    "# BLAST score distributions between Swiss-Prot variants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2609039a",
   "metadata": {},
   "outputs": [],
   "source": [
    "using Plots\n",
    "include(\"source/io.jl\")\n",
    "include(\"source/utils.jl\")\n",
    "include(\"source/derangement.jl\")\n",
    "include(\"source/alignment.jl\")\n",
    "include(\"source/blast.jl\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e22e726",
   "metadata": {},
   "source": [
    "### Shuffling and reversing \n",
    "The command `esl-shuffle input.fa > output.fa` shuffles each sequence in a FASTA record in `input.fa` and writes to `output.fa`. Using the `-r` flag, `esl-shuffle -r input.fa > output.fa` reverses each sequence. Generate the shuffled and reversed variants using EASEL.\n",
    "\n",
    "```\n",
    "esl-shuffle data/sprot.fa > data/shuf_sprot.fa\n",
    "esl-shuffle -r data/sprot.fa > data/rev_sprot.fa\n",
    "esl-shuffle -r data/shuf_sprot.fa > data/rev_shuf_sprot.fa\n",
    "\n",
    "esl-shuffle data/sprot_scop.fa > data/shuf_sprot_scop.fa\n",
    "esl-shuffle data/sprot_scop.fa -r > data/rev_sprot_scop.fa\n",
    "esl-shuffle -r data/shuf_sprot_scop.fa > data/rev_shuf_sprot_scop.fa\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75d5e661",
   "metadata": {},
   "source": [
    "Alternatively, generate the variants using Julia."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9edeb2c",
   "metadata": {},
   "outputs": [],
   "source": [
    "sprot = readsequences(\"data/sprot.fa\")\n",
    "rev_sprot = reverse.(sprot)\n",
    "writesequences(\"data/rev_sprot.fa\", rev_sprot)\n",
    "shuf_sprot = shufflefast.(sprot)\n",
    "writesequences(\"data/shuf_sprot.fa\", shuf_sprot)\n",
    "rev_shuf_sprot = reverse.(shuf_sprot)\n",
    "writesequences(\"data/rev_shuf_sprot.fa\", rev_shuf_sprot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c41758ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "sprot_scop = readsequences(\"data/sprot_scop.fa\")\n",
    "rev_sprot_scop = reverse.(sprot_scop)\n",
    "writesequences(\"data/rev_sprot_scop.fa\", rev_sprot_scop)\n",
    "shuf_sprot_scop = shufflefast.(sprot_scop)\n",
    "writesequences(\"data/shuf_sprot_scop.fa\", shuf_sprot_scop)\n",
    "rev_shuf_sprot_scop = reverse.(shuf_sprot_scop)\n",
    "writesequences(\"data/rev_shuf_sprot_scop.fa\", rev_shuf_sprot_scop)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b501dfac",
   "metadata": {},
   "source": [
    "### Classify Swiss-Prot sequences by their SCOP2 fold\n",
    "Download the SCOP2 classification table.\n",
    "\n",
    "```\n",
    "cd data\n",
    "wget http://scop.mrc-lmb.cam.ac.uk/files/scop-cla-latest.txt\n",
    "```\n",
    "\n",
    "If the cells for Figure 1 have not yet been run, run them now. The following cell requires `outputs/sprot_by_scop2sf.df`.\n",
    "\n",
    "Pair each sequence in `data/sprot_scop.fa` to its SCOP2 superfamily representative. Assign the Swiss-Prot sequence the fold of its SCOP2 superfamily representative and store these in a CSV."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "504edfdc",
   "metadata": {},
   "outputs": [],
   "source": [
    "scoptable = readscopclass(\"data/scop-cla-latest.txt\")\n",
    "foldlookup = Dict(Pair.(parse.(Int, scoptable.\"SF-DOMID\"), scoptable.\"CF\"));\n",
    "sprot_scop_registration = readframe(\"outputs/sprot_scop_registration.df\")\n",
    "sprot_scop_registration.\"fold\" = [foldlookup[scopid] for scopid=sprot_scop_registration.\"scop-id\"]\n",
    "writeframe(\"outputs/sprot_scop_registration.df\", sprot_scop_registration)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2436d8d",
   "metadata": {},
   "source": [
    "### Construct a derangement of the SCOP2 SF fold classification of Swiss-Prot.\n",
    "A derangement is a permutation that has no fixed points: **(1, 2, 3, 4) -> (3, 1, 4, 2)** is a derangement; **(1, 2, 3, 4) -> (3, 2, 4, 1)** is not. \n",
    "\n",
    "Construct a derangement $\\rho$ with respect to the fold classification such that all pairs of sequences $(s_i, s_{\\rho(i)})$ between Swiss-Prot and its reordering have distinct folds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75e3ac78",
   "metadata": {},
   "outputs": [],
   "source": [
    "sprot_scop_registration = readframe(\"outputs/sprot_scop_registration.df\")\n",
    "# the derangement is constructed from a vector of integers.\n",
    "# assign indices to each fold identified with a Swiss-Prot sequence. \n",
    "unique_folds = unique(sprot_scop_registration.\"fold\")\n",
    "indexlookup = Dict(Pair.(unique_folds, 1:length(unique_folds)))\n",
    "# then generate a vector of integers mapping sprot_scop index to fold index.\n",
    "sprot_scop_fold_indices = [indexlookup[fold] for fold=sprot_scop_registration.\"fold\"]\n",
    "# construct a derangenment of the sprot_scop fold sequence.\n",
    "ρ = derangement(sprot_scop_fold_indices, weighted=true)\n",
    "# and permute sprot by the same derangement.\n",
    "sprot = readfasta(\"data/sprot.fa\")\n",
    "sprot_scop = [sprot[i] for i=sprot_scop_registration.\"sprot-index\"]\n",
    "reorder_sprot_scop = sprot_scop[ρ]\n",
    "writefasta(\"data/reorder_sprot_scop.fa\", reorder_sprot_scop)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ecb60b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "### Paired alignment with BLAST\n",
    "With each variant constructed, run the full experiment structure using the Julia wrapper for `BLAST`.\n",
    "\n",
    "<img src=\"assets/diagram.png\" width=500 alt=\"The words sprot, shuf, rev, shufrev, and random are drawn in print. Four lines are drawn between these words, representing paired alignment between the sequence sets. A purple line connects sprot to random. A blue line connects sprot to rev. An orange line connects shuf to sprot. A red line connects shuf to shufrev.\"/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7422e252",
   "metadata": {},
   "outputs": [],
   "source": [
    "queries = [\n",
    "    \"data/sprot_scop.fa\", \n",
    "    \"data/sprot_scop.fa\", \n",
    "    \"data/sprot_scop.fa\", \n",
    "    \"data/shuf_sprot_scop.fa\"]\n",
    "references = [\n",
    "    \"data/reorder_sprot_scop.fa\", \n",
    "    \"data/shuf_sprot_scop.fa\", \n",
    "    \"data/rev_sprot_scop.fa\", \n",
    "    \"data/rev_shuf_sprot_scop.fa\"]\n",
    "name(x) = split(split(x, '/')[2], '.')[1]\n",
    "querynames = name.(queries)\n",
    "referencenames = name.(references)\n",
    "resultnames = [\"$(qname)↔$(rname)\" for (qname, rname)=zip(querynames, referencenames)]\n",
    "resultdirs = [\"data/$(name)\" for name=resultnames]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d16b556",
   "metadata": {},
   "outputs": [],
   "source": [
    "# this takes a while!\n",
    "search.([Pairwise()], resultdirs, queries, references; verbose=true, careful=true, ntasks=12)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2576dd7",
   "metadata": {},
   "source": [
    "## Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b82a95e",
   "metadata": {},
   "outputs": [],
   "source": [
    "using Plots\n",
    "using FASTX\n",
    "include(\"source/io.jl\")\n",
    "include(\"source/utils.jl\")\n",
    "function fig2plot!(x, label, color)\n",
    "    stephist!(x, label=label, color=color, bins=1:80)\n",
    "end\n",
    "colors = [\"purple\", \"orange\", \"blue\", \"red\"]\n",
    "results = readframe.(resultpaths)\n",
    "scores = [frame.score for frame=results]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46dc908e",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(xlabel=\"Score\",\n",
    "     ylabel=\"Alignments\",\n",
    "     legend=:outertop,\n",
    "     dpi=500)\n",
    "fig2plot!(scores[1], latexstring(BLAST(sprotscop, reorder(sprotscop))), colors[1])\n",
    "fig2plot!(scores[2], latexstring(BLAST(sprotscop, shuf(sprotscop))), colors[2])\n",
    "fig2plot!(scores[3], latexstring(BLAST(sprotscop, rev(sprotscop))), colors[3])\n",
    "fig2plot!(scores[4], latexstring(BLAST(shuf(sprotscop), rev(shuf(sprotscop)))), colors[4])\n",
    "savefig(\"figures/BLAST_SprotVariants.png\")\n",
    "plot!()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.8.5",
   "language": "julia",
   "name": "julia-1.8"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}