.template.cross-link_ms.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "%gencelloutputs\n",
    "Introduction to cross-linking {#cross-link_ms}\n",
    "============================="
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example we setup and evaluate a scoring function based on XL-MS data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "%%colabonly\n",
    "First, we need to install IMP:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%%colabonly\n",
    "!add-apt-repository -y ppa:salilab/ppa\n",
    "!apt install imp\n",
    "import sys, os, glob\n",
    "sys.path.append(os.path.dirname(glob.glob('/usr/lib/python*/dist-packages/IMP')[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "\n",
    "import IMP\n",
    "import IMP.pmi\n",
    "import IMP.pmi.topology\n",
    "import IMP.pmi.io\n",
    "import IMP.pmi.io.crosslink\n",
    "import IMP.pmi.restraints\n",
    "import IMP.pmi.restraints.crosslinking\n",
    "import ihm.cross_linkers\n",
    "\n",
    "m = IMP.Model()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we create the representation using PMI. We create two proteins:\n",
    "\n",
    " - ProtA with one bead, which coarse grains residues 1 to 10\n",
    " - ProtB with three beads, which coarse grains residues 1 to 10, 11 to 20 and 21 to 30."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "s = IMP.pmi.topology.System(m)\n",
    " \n",
    "st1 = s.create_state()\n",
    "protA = st1.create_molecule(\"ProtA\", \"G\" * 10, \"A\")\n",
    "protA.add_representation(resolutions=[10])\n",
    "protB = st1.create_molecule(\"ProtB\", \"G\" * 30, \"B\")\n",
    "protB.add_representation(resolutions=[10])\n",
    "hier = s.build()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To make it easier to see what's going on, we'll place the beads at fixed points in space:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "beads = IMP.atom.Selection(hier).get_selected_particles()\n",
    "print(beads)\n",
    "xyzs = [IMP.core.XYZ(b) for b in beads if IMP.core.XYZ.get_is_setup(b)]\n",
    "xyzs[0].set_coordinates((0,0,0))\n",
    "xyzs[1].set_coordinates((-40,0,0))\n",
    "xyzs[2].set_coordinates((0,0,0))\n",
    "xyzs[3].set_coordinates((40,0,0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we'll make some cross-links. The cross-link dataset is a comma separated value (CSV) file with at least the protein and the residue names (no spaces between commas):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "xldb='''Protein 1,Protein 2,Residue 1,Residue 2,UniqueID,Score\n",
    "ProtA,ProtB,1,10,1,1.0\n",
    "ProtA,ProtB,1,11,2,2.0\n",
    "ProtA,ProtB,1,21,3,2.0\n",
    "'''\n",
    "with open('xlinks.csv', 'w') as fh:\n",
    "    fh.write(xldb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we create a conversion map between internal keywords of cross-links features and the one in the file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "cldbkc = IMP.pmi.io.crosslink.CrossLinkDataBaseKeywordsConverter()\n",
    "cldbkc.set_protein1_key(\"Protein 1\")\n",
    "cldbkc.set_protein2_key(\"Protein 2\")\n",
    "cldbkc.set_residue1_key(\"Residue 1\")\n",
    "cldbkc.set_residue2_key(\"Residue 2\")\n",
    "\n",
    "# the unique_id_key and id_score_key are optional,\n",
    "# and they add features that will be explained below\n",
    "\n",
    "cldbkc.set_unique_id_key(\"UniqueID\")\n",
    "cldbkc.set_id_score_key(\"Score\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With this keyword interpreter, let's read the cross-link database:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "cldb = IMP.pmi.io.crosslink.CrossLinkDataBase(cldbkc)\n",
    "cldb.create_set_from_file(\"xlinks.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's check that the database looks ok:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "print(cldb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the database we can now set up the scoring function. Note the text generated. The program reports the nuisance particles associated to the cross-link (sigma and psi):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "xl = IMP.pmi.restraints.crosslinking.CrossLinkingMassSpectrometryRestraint(\n",
    "    root_hier=hier, database=cldb, length=21.0, slope=0.0,\n",
    "    resolution=1.0, label=\"XL\", linker=ihm.cross_linkers.dss)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can evaluate this restraint at the current system configuration:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "print(xl.rs.unprotected_evaluate(None))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot the score while moving ProtA bead wrt ProtB. First, we get the particle corresponding to ProtA:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "sel = IMP.atom.Selection(hier, molecule=\"ProtA\")\n",
    "pA, = sel.get_selected_particles()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can move ProtA on the x-axis:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "scores = []\n",
    "xs = []\n",
    "for i in range(-100, 100):\n",
    "    xs.append(float(i))\n",
    "    IMP.core.XYZ(pA).set_coordinates((i, 0, 0))\n",
    "    scores.append(xl.rs.unprotected_evaluate(None))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can make a simple plot of the scores against the x coordinate. This plot shows that the system has two minima:\n",
    " - one when ProtA is between ProtB:1-10 and ProtB:11-20;\n",
    " - and the other when ProtA is between ProtB:11-20 and ProtB:21-30"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import pylab\n",
    "\n",
    "pylab.plot(xs, scores)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plot is weird, so let's analyse what is going on.\n",
    "\n",
    "First let's simplify our dataset, by considering only the first cross-link. Let's filter by the UniqueID, creating a new database that contains only the second cross-link, namely UniqueID=2:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "from IMP.pmi.io.crosslink import FilterOperator\n",
    "import operator\n",
    "\n",
    "fo = FilterOperator(cldb.unique_id_key, operator.eq, \"2\")\n",
    "fcldb = cldb.filter(fo)\n",
    "\n",
    "print(fcldb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can create a new restraint based on this database and, as before, score while moving ProtA:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "xl1 = IMP.pmi.restraints.crosslinking.CrossLinkingMassSpectrometryRestraint(\n",
    "    root_hier=hier, database=fcldb, length=21.0, slope=0.0,\n",
    "    resolution=1.0, label=\"XL\", linker=ihm.cross_linkers.dss)\n",
    "\n",
    "scores=[]\n",
    "xs=[]\n",
    "for i in range(-100, 100):\n",
    "    xs.append(float(i))\n",
    "    IMP.core.XYZ(pA).set_coordinates((i, 0, 0))\n",
    "    scores.append(xl1.rs.unprotected_evaluate(None))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is clear that the restraint has a minimum when ProtA and ProtB:11-20 are close (namely when ProtA x is around 0). In fact, the restraint has a sigmoid shape:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "pylab.plot(xs, scores)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's play with the parameters sigma and psi to understand their roles. Let's get sigma first:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "sigma = xl1.sigma_dictionary[\"SIGMA\"][0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and let's vary its value between 1 and 20 to see what happens:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "scores_list = []\n",
    "xs_list = []\n",
    "for s in range(1, 20):\n",
    "    scores = []\n",
    "    xs = []\n",
    "    sigma.set_scale(float(s))\n",
    "    for i in range(-100, 100):\n",
    "        xs.append(float(i))\n",
    "        IMP.core.XYZ(pA).set_coordinates((i, 0, 0))\n",
    "        scores.append(xl1.rs.unprotected_evaluate(None))\n",
    "    scores_list.append(scores)\n",
    "    xs_list.append(xs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the plot, one can see sigma modulates both the slope of the sigmoid and the plateau of the minimum. This is because sigma is the structural uncertainty associated with the position of the cross-linked beads:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "for xs, scores in zip(xs_list, scores_list):\n",
    "    pylab.plot(xs, scores)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's get psi now (setting sigma back to 11), and vary its value between 0.01 and 0.5 to see what happens:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "sigma.set_scale(11)\n",
    "psi = xl1.psi_dictionary[\"PSI\"][0]\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "scores_list = []\n",
    "xs_list = []\n",
    "for s in np.linspace(0.01, 0.5, 10):\n",
    "    scores = []\n",
    "    xs = []\n",
    "    psi.set_scale(float(s))\n",
    "    for i in range(-100, 100):\n",
    "        xs.append(float(i))\n",
    "        IMP.core.XYZ(pA).set_coordinates((i, 0, 0))\n",
    "        scores.append(xl1.rs.unprotected_evaluate(None))\n",
    "    scores_list.append(scores)\n",
    "    xs_list.append(xs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting the values again, one can see psi modulates the plateau of the minimum and the maxima.\n",
    "This is because psi is the uncertainty associated with the cross-link observation:\n",
    "\n",
    " - When psi is 0.01, (low uncertainty) there is a big score difference between a satisfied cross-link (x=0) and a violated one (x=100, x=-100)\n",
    " - When psi is 0.5, (high uncertainty) the score is flat, and there is no difference between violated and satisfied cross-links"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "for xs, scores in zip(xs_list, scores_list):\n",
    "    pylab.plot(xs, scores)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's move the parameters when we consider all three cross-links. First, sigma:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "sigma = xl.sigma_dictionary[\"SIGMA\"][0]\n",
    "\n",
    "scores_list = []\n",
    "xs_list = []\n",
    "for s in range(1, 20):\n",
    "    scores = []\n",
    "    xs = []\n",
    "    sigma.set_scale(float(s))\n",
    "    for i in range(-100, 100):\n",
    "        xs.append(float(i))\n",
    "        IMP.core.XYZ(pA).set_coordinates((i, 0, 0))\n",
    "        scores.append(xl.rs.unprotected_evaluate(None))\n",
    "    scores_list.append(scores)\n",
    "    xs_list.append(xs)\n",
    "\n",
    "for xs, scores in zip(xs_list, scores_list):\n",
    "    pylab.plot(xs, scores)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, for psi:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "sigma.set_scale(11)\n",
    "psi = xl.psi_dictionary[\"PSI\"][0]\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "scores_list = []\n",
    "xs_list = []\n",
    "for s in np.linspace(0.01, 0.5, 10):\n",
    "    scores = []\n",
    "    xs = []\n",
    "    psi.set_scale(float(s))\n",
    "    for i in range(-100, 100):\n",
    "        xs.append(float(i))\n",
    "        IMP.core.XYZ(pA).set_coordinates((i, 0, 0))\n",
    "        scores.append(xl.rs.unprotected_evaluate(None))\n",
    "    scores_list.append(scores)\n",
    "    xs_list.append(xs)\n",
    "\n",
    "for xs, scores in zip(xs_list, scores_list):\n",
    "    pylab.plot(xs, scores)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can try to optimize the values of PSI and SIGMA, and see what is the best scoring value, fixing the coordinate of ProtA to a minimum:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "import math\n",
    "import numpy as np\n",
    "\n",
    "sel = IMP.atom.Selection(hier, molecule=\"ProtB\")\n",
    "pB1, pB2, pB3 = sel.get_selected_particles()\n",
    "\n",
    "IMP.core.XYZ(pA).set_coordinates((0,0,0))\n",
    "IMP.core.XYZ(pB1).set_coordinates((0,20,0))\n",
    "IMP.core.XYZ(pB2).set_coordinates((-20*math.sqrt(3)/2,-20/2,0))\n",
    "IMP.core.XYZ(pB3).set_coordinates((20*math.sqrt(3)/2,-20/2,0))\n",
    "\n",
    "scores = []\n",
    "psis = []\n",
    "sigmas = []\n",
    "for p in np.linspace(0.01, 0.5, 100):\n",
    "    psi.set_scale(p)\n",
    "    for s in np.linspace(1, 40, 50):\n",
    "        psis.append(p)\n",
    "        sigmas.append(s)\n",
    "        sigma.set_scale(s)\n",
    "        scores.append(xl.rs.unprotected_evaluate(None))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is a minimum when PSI is close to zero and sigma is between 0 and 10:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.scatter(psis, sigmas, c=scores, s=30, edgecolor=[])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "Let's move ProtA away, so that any cross-link is satisfied:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "IMP.core.XYZ(pA).set_coordinates((100, 100, 100))\n",
    "scores = []\n",
    "psis = []\n",
    "sigmas = []\n",
    "for p in np.linspace(0.01, 0.5, 100):\n",
    "    psi.set_scale(p)\n",
    "    for s in np.linspace(1, 40, 50):\n",
    "        psis.append(p)\n",
    "        sigmas.append(s)\n",
    "        sigma.set_scale(s)\n",
    "        scores.append(xl.rs.unprotected_evaluate(None))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The minimum is at Psi=0.5, irrespective of the value of Sigma:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.scatter(psis, sigmas, c=scores, s=30, edgecolor=[])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Further reading\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "%%colabexclude\n",
    "Ambiguity can also be an issue when dealing with cross-links. See [this followup tutorial](@file cross-link_ms-ambiguity) for further information."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "%%colabonly\n",
    "Ambiguity can also be an issue when dealing with cross-links. See [this followup tutorial](@file cross-link_ms-ambiguity-colab) for further information."
   ]
  }
 ],
 "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}