Created
August 8, 2022 14:58
-
-
Save alisterburt/706cbdff00e7011064dbf09589bbe972 to your computer and use it in GitHub Desktop.
ribosome tilt-series simulation
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| { | |
| "cells": [ | |
| { | |
| "cell_type": "code", | |
| "execution_count": 1, | |
| "metadata": { | |
| "collapsed": true | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "from pathlib import Path\n", | |
| "\n", | |
| "import numpy as np\n", | |
| "import pandas as pd\n", | |
| "import gemmi\n", | |
| "import mrcfile\n", | |
| "import starfile\n", | |
| "import einops\n", | |
| "import eulerangles\n", | |
| "from scipy.stats import special_ortho_group\n", | |
| "from scipy.spatial.transform import Rotation as R\n" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 2, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": "\u001B[1m<\u001B[0m\u001B[1;95mgemmi.Model\u001B[0m\u001B[39m \u001B[0m\u001B[1;36m1\u001B[0m\u001B[39m with \u001B[0m\u001B[1;36m89\u001B[0m\u001B[39m \u001B[0m\u001B[1;35mchain\u001B[0m\u001B[1;39m(\u001B[0m\u001B[39ms\u001B[0m\u001B[1;39m)\u001B[0m\u001B[1m>\u001B[0m\n", | |
| "text/html": "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"><span style=\"font-weight: bold\"><</span><span style=\"color: #ff00ff; text-decoration-color: #ff00ff; font-weight: bold\">gemmi.Model</span><span style=\"color: #000000; text-decoration-color: #000000\"> </span><span style=\"color: #008080; text-decoration-color: #008080; font-weight: bold\">1</span><span style=\"color: #000000; text-decoration-color: #000000\"> with </span><span style=\"color: #008080; text-decoration-color: #008080; font-weight: bold\">89</span><span style=\"color: #000000; text-decoration-color: #000000\"> </span><span style=\"color: #800080; text-decoration-color: #800080; font-weight: bold\">chain</span><span style=\"color: #000000; text-decoration-color: #000000; font-weight: bold\">(</span><span style=\"color: #000000; text-decoration-color: #000000\">s</span><span style=\"color: #000000; text-decoration-color: #000000; font-weight: bold\">)</span><span style=\"font-weight: bold\">></span>\n</pre>\n" | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "pdb = '4v6x-ribo.cif'\n", | |
| "structure = gemmi.read_structure(pdb)\n", | |
| "model = structure[0]\n", | |
| "model" | |
| ], | |
| "metadata": { | |
| "collapsed": false, | |
| "pycharm": { | |
| "name": "#%%\n" | |
| } | |
| } | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 3, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": "\u001B[1m(\u001B[0m\u001B[1;36m13338\u001B[0m, \u001B[1;36m3\u001B[0m\u001B[1m)\u001B[0m\n", | |
| "text/html": "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"><span style=\"font-weight: bold\">(</span><span style=\"color: #008080; text-decoration-color: #008080; font-weight: bold\">13338</span>, <span style=\"color: #008080; text-decoration-color: #008080; font-weight: bold\">3</span><span style=\"font-weight: bold\">)</span>\n</pre>\n" | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "ca_coords = np.array([\n", | |
| " [cra.atom.pos.x, cra.atom.pos.y, cra.atom.pos.z]\n", | |
| " for cra in model.all()\n", | |
| " if cra.atom.name == 'CA'\n", | |
| "])\n", | |
| "ca_coords.shape" | |
| ], | |
| "metadata": { | |
| "collapsed": false, | |
| "pycharm": { | |
| "name": "#%%\n" | |
| } | |
| } | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 4, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": "\u001B[1m(\u001B[0m\u001B[1;36m13338\u001B[0m, \u001B[1;36m3\u001B[0m\u001B[1m)\u001B[0m\n", | |
| "text/html": "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"><span style=\"font-weight: bold\">(</span><span style=\"color: #008080; text-decoration-color: #008080; font-weight: bold\">13338</span>, <span style=\"color: #008080; text-decoration-color: #008080; font-weight: bold\">3</span><span style=\"font-weight: bold\">)</span>\n</pre>\n" | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "ca_coords_10apx = ca_coords / 10\n", | |
| "ca_coords_10apx_centered = ca_coords_10apx - np.mean(ca_coords_10apx, axis=0)\n", | |
| "assert np.allclose(np.mean(ca_coords_10apx_centered, axis=0), [0, 0, 0])\n", | |
| "ca_coords_10apx_centered.shape" | |
| ], | |
| "metadata": { | |
| "collapsed": false, | |
| "pycharm": { | |
| "name": "#%%\n" | |
| } | |
| } | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 5, | |
| "outputs": [], | |
| "source": [ | |
| "model_coords = ca_coords_10apx_centered" | |
| ], | |
| "metadata": { | |
| "collapsed": false, | |
| "pycharm": { | |
| "name": "#%%\n" | |
| } | |
| } | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 6, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": "\u001B[1m(\u001B[0m\u001B[1;35marray\u001B[0m\u001B[1m(\u001B[0m\u001B[1m[\u001B[0m\u001B[1;36m-13.74219601\u001B[0m, \u001B[1;36m-12.59160924\u001B[0m, \u001B[1;36m-17.41221623\u001B[0m\u001B[1m]\u001B[0m\u001B[1m)\u001B[0m, \u001B[1;35marray\u001B[0m\u001B[1m(\u001B[0m\u001B[1m[\u001B[0m\u001B[1;36m12.15490399\u001B[0m, \u001B[1;36m14.24029076\u001B[0m, \u001B[1;36m12.01938377\u001B[0m\u001B[1m]\u001B[0m\u001B[1m)\u001B[0m\u001B[1m)\u001B[0m\n", | |
| "text/html": "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"><span style=\"font-weight: bold\">(</span><span style=\"color: #800080; text-decoration-color: #800080; font-weight: bold\">array</span><span style=\"font-weight: bold\">([</span><span style=\"color: #008080; text-decoration-color: #008080; font-weight: bold\">-13.74219601</span>, <span style=\"color: #008080; text-decoration-color: #008080; font-weight: bold\">-12.59160924</span>, <span style=\"color: #008080; text-decoration-color: #008080; font-weight: bold\">-17.41221623</span><span style=\"font-weight: bold\">])</span>, <span style=\"color: #800080; text-decoration-color: #800080; font-weight: bold\">array</span><span style=\"font-weight: bold\">([</span><span style=\"color: #008080; text-decoration-color: #008080; font-weight: bold\">12.15490399</span>, <span style=\"color: #008080; text-decoration-color: #008080; font-weight: bold\">14.24029076</span>, <span style=\"color: #008080; text-decoration-color: #008080; font-weight: bold\">12.01938377</span><span style=\"font-weight: bold\">]))</span>\n</pre>\n" | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "np.min(model_coords, axis=0), np.max(model_coords, axis=0)" | |
| ], | |
| "metadata": { | |
| "collapsed": false, | |
| "pycharm": { | |
| "name": "#%%\n" | |
| } | |
| } | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 7, | |
| "outputs": [], | |
| "source": [ | |
| "# generate 20 ribos in a 500^3 box centered on 0,0,0" | |
| ], | |
| "metadata": { | |
| "collapsed": false, | |
| "pycharm": { | |
| "name": "#%%\n" | |
| } | |
| } | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 8, | |
| "outputs": [], | |
| "source": [ | |
| "n = 20\n", | |
| "model_positions = np.random.uniform(low=-200, high=200, size=(n, 3))\n", | |
| "model_orientations = special_ortho_group.rvs(dim=3, size=n)" | |
| ], | |
| "metadata": { | |
| "collapsed": false, | |
| "pycharm": { | |
| "name": "#%%\n" | |
| } | |
| } | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 9, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": "\u001B[1m(\u001B[0m\u001B[1;36m20\u001B[0m, \u001B[1;36m3\u001B[0m, \u001B[1;36m3\u001B[0m\u001B[1m)\u001B[0m\n", | |
| "text/html": "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"><span style=\"font-weight: bold\">(</span><span style=\"color: #008080; text-decoration-color: #008080; font-weight: bold\">20</span>, <span style=\"color: #008080; text-decoration-color: #008080; font-weight: bold\">3</span>, <span style=\"color: #008080; text-decoration-color: #008080; font-weight: bold\">3</span><span style=\"font-weight: bold\">)</span>\n</pre>\n" | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "model_orientations.shape\n" | |
| ], | |
| "metadata": { | |
| "collapsed": false, | |
| "pycharm": { | |
| "name": "#%%\n" | |
| } | |
| } | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 10, | |
| "outputs": [], | |
| "source": [ | |
| "model_orientations = model_orientations.reshape((-1, 1, 3, 3))" | |
| ], | |
| "metadata": { | |
| "collapsed": false, | |
| "pycharm": { | |
| "name": "#%%\n" | |
| } | |
| } | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 11, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": "\u001B[1m(\u001B[0m\u001B[1;36m13338\u001B[0m, \u001B[1;36m20\u001B[0m, \u001B[1;36m3\u001B[0m\u001B[1m)\u001B[0m\n", | |
| "text/html": "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"><span style=\"font-weight: bold\">(</span><span style=\"color: #008080; text-decoration-color: #008080; font-weight: bold\">13338</span>, <span style=\"color: #008080; text-decoration-color: #008080; font-weight: bold\">20</span>, <span style=\"color: #008080; text-decoration-color: #008080; font-weight: bold\">3</span><span style=\"font-weight: bold\">)</span>\n</pre>\n" | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "oriented_model_coordinates = model_orientations @ model_coords.reshape((-1, 3, 1))\n", | |
| "oriented_model_coordinates = np.squeeze(oriented_model_coordinates).transpose(1, 0, 2)\n", | |
| "oriented_model_coordinates.shape" | |
| ], | |
| "metadata": { | |
| "collapsed": false, | |
| "pycharm": { | |
| "name": "#%%\n" | |
| } | |
| } | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 12, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": "\u001B[1m(\u001B[0m\u001B[1m(\u001B[0m\u001B[1;36m20\u001B[0m, \u001B[1;36m3\u001B[0m\u001B[1m)\u001B[0m, \u001B[1m(\u001B[0m\u001B[1;36m13338\u001B[0m, \u001B[1;36m20\u001B[0m, \u001B[1;36m3\u001B[0m\u001B[1m)\u001B[0m\u001B[1m)\u001B[0m\n", | |
| "text/html": "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"><span style=\"font-weight: bold\">((</span><span style=\"color: #008080; text-decoration-color: #008080; font-weight: bold\">20</span>, <span style=\"color: #008080; text-decoration-color: #008080; font-weight: bold\">3</span><span style=\"font-weight: bold\">)</span>, <span style=\"font-weight: bold\">(</span><span style=\"color: #008080; text-decoration-color: #008080; font-weight: bold\">13338</span>, <span style=\"color: #008080; text-decoration-color: #008080; font-weight: bold\">20</span>, <span style=\"color: #008080; text-decoration-color: #008080; font-weight: bold\">3</span><span style=\"font-weight: bold\">))</span>\n</pre>\n" | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "positioned_oriented_model_coordinates = oriented_model_coordinates + model_positions\n", | |
| "model_positions.shape, positioned_oriented_model_coordinates.shape" | |
| ], | |
| "metadata": { | |
| "collapsed": false, | |
| "pycharm": { | |
| "name": "#%%\n" | |
| } | |
| } | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 13, | |
| "outputs": [], | |
| "source": [ | |
| "final_coordinates = positioned_oriented_model_coordinates.reshape((-1, 3))" | |
| ], | |
| "metadata": { | |
| "collapsed": false, | |
| "pycharm": { | |
| "name": "#%%\n" | |
| } | |
| } | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 14, | |
| "outputs": [], | |
| "source": [ | |
| "def render_on_grid_3d(coords: np.ndarray) -> np.ndarray:\n", | |
| " b = np.linspace(-250, 250, num=501, endpoint=True)\n", | |
| " image, _ = np.histogramdd(coords, bins=[b, b, b])\n", | |
| " return image\n", | |
| "\n", | |
| "\n", | |
| "def render_on_grid_2d(coords: np.ndarray) -> np.ndarray:\n", | |
| " b = np.linspace(-250.5, 250.5, num=501, endpoint=True)\n", | |
| " image, _ = np.histogramdd(coords[..., :2], bins=[b, b])\n", | |
| " return image" | |
| ], | |
| "metadata": { | |
| "collapsed": false, | |
| "pycharm": { | |
| "name": "#%%\n" | |
| } | |
| } | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 15, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": "\u001B[1m(\u001B[0m\u001B[1;36m500\u001B[0m, \u001B[1;36m500\u001B[0m\u001B[1m)\u001B[0m\n", | |
| "text/html": "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"><span style=\"font-weight: bold\">(</span><span style=\"color: #008080; text-decoration-color: #008080; font-weight: bold\">500</span>, <span style=\"color: #008080; text-decoration-color: #008080; font-weight: bold\">500</span><span style=\"font-weight: bold\">)</span>\n</pre>\n" | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "test = render_on_grid_2d(final_coordinates)\n", | |
| "test.shape" | |
| ], | |
| "metadata": { | |
| "collapsed": false, | |
| "pycharm": { | |
| "name": "#%%\n" | |
| } | |
| } | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 16, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": "\u001B[1m<\u001B[0m\u001B[1;95mImage\u001B[0m\u001B[39m layer \u001B[0m\u001B[32m'test'\u001B[0m\u001B[39m at \u001B[0m\u001B[1;36m0x7ffdf2fd7c10\u001B[0m\u001B[1m>\u001B[0m\n", | |
| "text/html": "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"><span style=\"font-weight: bold\"><</span><span style=\"color: #ff00ff; text-decoration-color: #ff00ff; font-weight: bold\">Image</span><span style=\"color: #000000; text-decoration-color: #000000\"> layer </span><span style=\"color: #008000; text-decoration-color: #008000\">'test'</span><span style=\"color: #000000; text-decoration-color: #000000\"> at </span><span style=\"color: #008080; text-decoration-color: #008080; font-weight: bold\">0x7ffdf2fd7c10</span><span style=\"font-weight: bold\">></span>\n</pre>\n" | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "import napari\n", | |
| "\n", | |
| "viewer = napari.Viewer()\n", | |
| "viewer.add_image(test)" | |
| ], | |
| "metadata": { | |
| "collapsed": false, | |
| "pycharm": { | |
| "name": "#%%\n" | |
| } | |
| } | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 17, | |
| "outputs": [], | |
| "source": [ | |
| "def Rx(angles_degrees: np.ndarray) -> np.ndarray:\n", | |
| " \"\"\"Affine matrix for a rotation around the X-axis.\"\"\"\n", | |
| " angles_degrees = np.asarray(angles_degrees).reshape(-1)\n", | |
| " c = np.cos(np.deg2rad(angles_degrees))\n", | |
| " s = np.sin(np.deg2rad(angles_degrees))\n", | |
| " matrices = einops.repeat(\n", | |
| " np.eye(4), 'i j -> n i j', n=len(angles_degrees)\n", | |
| " )\n", | |
| " matrices[:, 1, 1] = c\n", | |
| " matrices[:, 1, 2] = -s\n", | |
| " matrices[:, 2, 1] = s\n", | |
| " matrices[:, 2, 2] = c\n", | |
| " return np.squeeze(matrices)\n", | |
| "\n", | |
| "\n", | |
| "def Ry(angles_degrees: np.ndarray) -> np.ndarray:\n", | |
| " \"\"\"Affine matrix for a rotation around the Y-axis.\"\"\"\n", | |
| " angles_degrees = np.asarray(angles_degrees).reshape(-1)\n", | |
| " c = np.cos(np.deg2rad(angles_degrees))\n", | |
| " s = np.sin(np.deg2rad(angles_degrees))\n", | |
| " matrices = einops.repeat(\n", | |
| " np.eye(4), 'i j -> n i j', n=len(angles_degrees)\n", | |
| " )\n", | |
| " matrices[:, 0, 0] = c\n", | |
| " matrices[:, 0, 2] = s\n", | |
| " matrices[:, 2, 0] = -s\n", | |
| " matrices[:, 2, 2] = c\n", | |
| " return np.squeeze(matrices)\n", | |
| "\n", | |
| "\n", | |
| "def Rz(angles_degrees: float) -> np.ndarray:\n", | |
| " \"\"\"Affine matrix for a rotation around the Z-axis.\"\"\"\n", | |
| " angle_degrees = np.asarray(angles_degrees).reshape(-1)\n", | |
| " c = np.cos(np.deg2rad(angle_degrees))\n", | |
| " s = np.sin(np.deg2rad(angle_degrees))\n", | |
| " matrices = einops.repeat(\n", | |
| " np.eye(4), 'i j -> n i j', n=len(angle_degrees)\n", | |
| " )\n", | |
| " matrices[:, 0, 0] = c\n", | |
| " matrices[:, 0, 1] = -s\n", | |
| " matrices[:, 1, 0] = s\n", | |
| " matrices[:, 1, 1] = c\n", | |
| " return np.squeeze(matrices)\n", | |
| "\n", | |
| "\n", | |
| "def S(shifts: np.ndarray) -> np.ndarray:\n", | |
| " \"\"\"Affine matrices for shifts.\n", | |
| " Shifts supplied can be 2D or 3D.\n", | |
| " \"\"\"\n", | |
| " shifts = np.asarray(shifts, dtype=float)\n", | |
| " if shifts.shape[-1] == 2:\n", | |
| " shifts = _promote_2d_to_3d(shifts)\n", | |
| " shifts = np.array(shifts).reshape((-1, 3))\n", | |
| " matrices = einops.repeat(np.eye(4), 'i j -> n i j', n=shifts.shape[0])\n", | |
| " matrices[:, 0:3, 3] = shifts\n", | |
| " return np.squeeze(matrices)\n", | |
| "\n", | |
| "\n", | |
| "def _promote_2d_to_3d(shifts: np.ndarray) -> np.ndarray:\n", | |
| " \"\"\"Promote 2D vectors to 3D with zeros in the last dimension.\"\"\"\n", | |
| " shifts = np.asarray(shifts).reshape(-1, 2)\n", | |
| " shifts = np.c_[shifts, np.zeros(shifts.shape[0])]\n", | |
| " return np.squeeze(shifts)\n", | |
| "\n", | |
| "def homogenise_coordinates(coords: np.ndarray) -> np.ndarray:\n", | |
| " return np.c_[coords, np.ones(len(coords))]\n" | |
| ], | |
| "metadata": { | |
| "collapsed": false, | |
| "pycharm": { | |
| "name": "#%%\n" | |
| } | |
| } | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 18, | |
| "outputs": [], | |
| "source": [ | |
| "# simulate tilt-series\n", | |
| "def simulate_tilt_image(coords, rx=0, ry=0, rz=0, dx=0, dy=0):\n", | |
| " coords = homogenise_coordinates(coords)\n", | |
| " transformation = S([dx, dy, 0]) @ Rz(rz) @ Ry(ry) @ Rx(rx)\n", | |
| " transformed_coords = np.squeeze(transformation @ coords.reshape((-1, 4, 1)))\n", | |
| " return render_on_grid_2d(transformed_coords)" | |
| ], | |
| "metadata": { | |
| "collapsed": false, | |
| "pycharm": { | |
| "name": "#%%\n" | |
| } | |
| } | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 19, | |
| "outputs": [], | |
| "source": [ | |
| "test_0 = simulate_tilt_image(final_coordinates, ry=0)\n", | |
| "test_60 = simulate_tilt_image(final_coordinates, ry=5)" | |
| ], | |
| "metadata": { | |
| "collapsed": false, | |
| "pycharm": { | |
| "name": "#%%\n" | |
| } | |
| } | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 20, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": "\u001B[1m<\u001B[0m\u001B[1;95mImage\u001B[0m\u001B[39m layer \u001B[0m\u001B[32m'test_60'\u001B[0m\u001B[39m at \u001B[0m\u001B[1;36m0x7ffdd16103a0\u001B[0m\u001B[1m>\u001B[0m\n", | |
| "text/html": "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"><span style=\"font-weight: bold\"><</span><span style=\"color: #ff00ff; text-decoration-color: #ff00ff; font-weight: bold\">Image</span><span style=\"color: #000000; text-decoration-color: #000000\"> layer </span><span style=\"color: #008000; text-decoration-color: #008000\">'test_60'</span><span style=\"color: #000000; text-decoration-color: #000000\"> at </span><span style=\"color: #008080; text-decoration-color: #008080; font-weight: bold\">0x7ffdd16103a0</span><span style=\"font-weight: bold\">></span>\n</pre>\n" | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "viewer = napari.Viewer()\n", | |
| "viewer.add_image(test_0)\n", | |
| "viewer.add_image(test_60)" | |
| ], | |
| "metadata": { | |
| "collapsed": false, | |
| "pycharm": { | |
| "name": "#%%\n" | |
| } | |
| } | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 21, | |
| "outputs": [], | |
| "source": [ | |
| "tilt_series = np.stack(\n", | |
| " [\n", | |
| " simulate_tilt_image(final_coordinates, ry=r)\n", | |
| " for r\n", | |
| " in np.arange(-90, 90, 3)\n", | |
| " ]\n", | |
| ")" | |
| ], | |
| "metadata": { | |
| "collapsed": false, | |
| "pycharm": { | |
| "name": "#%%\n" | |
| } | |
| } | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 22, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": "\u001B[1;36m60\u001B[0m\n", | |
| "text/html": "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"><span style=\"color: #008080; text-decoration-color: #008080; font-weight: bold\">60</span>\n</pre>\n" | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "len(np.arange(-90, 90, 3))" | |
| ], | |
| "metadata": { | |
| "collapsed": false, | |
| "pycharm": { | |
| "name": "#%%\n" | |
| } | |
| } | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 23, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": "\u001B[1m(\u001B[0m\u001B[1;36m20\u001B[0m, \u001B[1;36m60\u001B[0m, \u001B[1;36m2\u001B[0m\u001B[1m)\u001B[0m\n", | |
| "text/html": "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"><span style=\"font-weight: bold\">(</span><span style=\"color: #008080; text-decoration-color: #008080; font-weight: bold\">20</span>, <span style=\"color: #008080; text-decoration-color: #008080; font-weight: bold\">60</span>, <span style=\"color: #008080; text-decoration-color: #008080; font-weight: bold\">2</span><span style=\"font-weight: bold\">)</span>\n</pre>\n" | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "tilt_series = tilt_series.transpose(0, 2, 1)\n", | |
| "particle_positions = model_positions + np.array([250, 250, 250])\n", | |
| "\n", | |
| "transformations = S([250, 250]) @ Ry(np.arange(-90, 90, 3)) @ S([-250, -250, -250])\n", | |
| "projected_positions = np.squeeze(transformations @ homogenise_coordinates(particle_positions).reshape((-1, 1, 4, 1)))[..., :2]\n", | |
| "projected_positions.shape" | |
| ], | |
| "metadata": { | |
| "collapsed": false, | |
| "pycharm": { | |
| "name": "#%%\n" | |
| } | |
| } | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 24, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": "\u001B[1m(\u001B[0m\u001B[1;36m1200\u001B[0m, \u001B[1;36m3\u001B[0m\u001B[1m)\u001B[0m\n", | |
| "text/html": "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"><span style=\"font-weight: bold\">(</span><span style=\"color: #008080; text-decoration-color: #008080; font-weight: bold\">1200</span>, <span style=\"color: #008080; text-decoration-color: #008080; font-weight: bold\">3</span><span style=\"font-weight: bold\">)</span>\n</pre>\n" | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "visualised_positions = np.zeros((n, 60, 3))\n", | |
| "visualised_positions[..., 0] = np.arange(60)\n", | |
| "visualised_positions[..., 1:] = projected_positions[..., ::-1]\n", | |
| "visualised_positions = visualised_positions.reshape((-1, 3))\n", | |
| "visualised_positions.shape" | |
| ], | |
| "metadata": { | |
| "collapsed": false, | |
| "pycharm": { | |
| "name": "#%%\n" | |
| } | |
| } | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 25, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": "\u001B[1m<\u001B[0m\u001B[1;95mPoints\u001B[0m\u001B[39m layer \u001B[0m\u001B[32m'visualised_positions'\u001B[0m\u001B[39m at \u001B[0m\u001B[1;36m0x7ffdc368b370\u001B[0m\u001B[1m>\u001B[0m\n", | |
| "text/html": "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"><span style=\"font-weight: bold\"><</span><span style=\"color: #ff00ff; text-decoration-color: #ff00ff; font-weight: bold\">Points</span><span style=\"color: #000000; text-decoration-color: #000000\"> layer </span><span style=\"color: #008000; text-decoration-color: #008000\">'visualised_positions'</span><span style=\"color: #000000; text-decoration-color: #000000\"> at </span><span style=\"color: #008080; text-decoration-color: #008080; font-weight: bold\">0x7ffdc368b370</span><span style=\"font-weight: bold\">></span>\n</pre>\n" | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "viewer = napari.Viewer()\n", | |
| "viewer.add_image(tilt_series)\n", | |
| "viewer.add_points(visualised_positions)" | |
| ], | |
| "metadata": { | |
| "collapsed": false, | |
| "pycharm": { | |
| "name": "#%%\n" | |
| } | |
| } | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 26, | |
| "outputs": [], | |
| "source": [ | |
| "output_dir = Path('tilt_images')\n", | |
| "output_dir.mkdir(exist_ok=True)\n", | |
| "for image, tilt_angle in zip(tilt_series, np.arange(-90, 90, 3)):\n", | |
| " mrcfile.write(output_dir / f'simulated_tilt_image_{tilt_angle:.2f}.mrc', tilt_series.astype(np.float32), voxel_size=10, overwrite=True)" | |
| ], | |
| "metadata": { | |
| "collapsed": false, | |
| "pycharm": { | |
| "name": "#%%\n" | |
| } | |
| } | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 27, | |
| "outputs": [], | |
| "source": [ | |
| "global_star_file = 'aligned_tilt_series.star'\n", | |
| "individual_star_file = 'alignments_simulated_ts.star'\n", | |
| "particle_star_file = 'particles.star'" | |
| ], | |
| "metadata": { | |
| "collapsed": false, | |
| "pycharm": { | |
| "name": "#%%\n" | |
| } | |
| } | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 28, | |
| "outputs": [], | |
| "source": [ | |
| "global_data = pd.DataFrame({\n", | |
| " 'rlnTomoName': ['TS_01'],\n", | |
| " 'rlnTomoTiltSeriesStarFile': [individual_star_file],\n", | |
| " 'rlnVoltage': [300],\n", | |
| " 'rlnSphericalAberration': [2.7],\n", | |
| " 'rlnAmplitudeContrast': [0.1],\n", | |
| " 'rlnMicrographOriginalPixelSize': [10],\n", | |
| " 'rlnTomoHand': [1],\n", | |
| " 'rlnTomoTiltSeriesPixelSize': [10]\n", | |
| "})\n", | |
| "\n", | |
| "starfile.write({'global': global_data}, global_star_file, overwrite=True)" | |
| ], | |
| "metadata": { | |
| "collapsed": false, | |
| "pycharm": { | |
| "name": "#%%\n" | |
| } | |
| } | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 29, | |
| "outputs": [], | |
| "source": [ | |
| "alignment_data = pd.DataFrame({\n", | |
| " 'rlnMicrographName': [output_dir / f'simulated_tilt_image_{tilt_angle:.2f}.mrc' for tilt_angle in np.arange(-90, 90, 3)],\n", | |
| " 'rlnMicrographPreExposure': [0 for _ in range(len(np.arange(-90, 90, 3)))],\n", | |
| " 'rlnDefocusU': [0 for _ in range(len(np.arange(-90, 90, 3)))],\n", | |
| " 'rlnDefocusV': [0 for _ in range(len(np.arange(-90, 90, 3)))],\n", | |
| " 'rlnDefocusAngle': [0 for _ in range(len(np.arange(-90, 90, 3)))],\n", | |
| "})\n", | |
| "\n", | |
| "projection_matrix_labels = [f'rlnTomoProj{ax}' for ax in 'XYZW']\n", | |
| "for idx, label in enumerate(projection_matrix_labels):\n", | |
| " rows = transformations[:, idx, :]\n", | |
| " alignment_data[label] = [\n", | |
| " f'[{r[0]:.13g},{r[1]:.13g},{r[2]:.13g},{r[3]:.13g}]'\n", | |
| " for r in rows\n", | |
| " ]\n", | |
| "\n", | |
| "starfile.write({'TS_01': alignment_data}, individual_star_file, overwrite=True)" | |
| ], | |
| "metadata": { | |
| "collapsed": false, | |
| "pycharm": { | |
| "name": "#%%\n" | |
| } | |
| } | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 30, | |
| "outputs": [], | |
| "source": [ | |
| "particle_optics = pd.DataFrame({\n", | |
| " 'rlnOpticsGroup': [1],\n", | |
| " 'rlnOpticsGroupName': ['optics'],\n", | |
| " 'rlnSphericalAberration': [2.7],\n", | |
| " 'rlnVoltage': [300],\n", | |
| " 'rlnTomoTiltSeriesPixelSize': [10],\n", | |
| "})\n", | |
| "\n", | |
| "particle_eulers = eulerangles.matrix2euler(\n", | |
| " rotation_matrices=np.linalg.pinv(model_orientations),\n", | |
| " axes='zyz',\n", | |
| " intrinsic=True,\n", | |
| " right_handed_rotation=True,\n", | |
| " )\n", | |
| "particle_data = pd.DataFrame({\n", | |
| " 'rlnTomoName': ['TS_01' for _ in range(len(particle_positions))],\n", | |
| " 'rlnTomoParticleId': [i for i in range(len(particle_positions))],\n", | |
| " 'rlnCoordinateX': particle_positions[:, 0],\n", | |
| " 'rlnCoordinateY': particle_positions[:, 1],\n", | |
| " 'rlnCoordinateZ': particle_positions[:, 2],\n", | |
| " 'rlnAngleRot': particle_eulers[:, 0],\n", | |
| " 'rlnAngleTilt': particle_eulers[:, 1],\n", | |
| " 'rlnAnglePsi': particle_eulers[:, 2],\n", | |
| "})\n", | |
| "\n", | |
| "starfile.write({'optics': particle_optics, 'particles': particle_data}, particle_star_file, overwrite=True)" | |
| ], | |
| "metadata": { | |
| "collapsed": false, | |
| "pycharm": { | |
| "name": "#%%\n" | |
| } | |
| } | |
| } | |
| ], | |
| "metadata": { | |
| "kernelspec": { | |
| "display_name": "Python 3", | |
| "language": "python", | |
| "name": "python3" | |
| }, | |
| "language_info": { | |
| "codemirror_mode": { | |
| "name": "ipython", | |
| "version": 2 | |
| }, | |
| "file_extension": ".py", | |
| "mimetype": "text/x-python", | |
| "name": "python", | |
| "nbconvert_exporter": "python", | |
| "pygments_lexer": "ipython2", | |
| "version": "2.7.6" | |
| } | |
| }, | |
| "nbformat": 4, | |
| "nbformat_minor": 0 | |
| } |
Author
Only because they are missing CTF info? No CTF is applied to the simulation so it doesn’t make sense to add it
this was written to test some new features :)
for rln40
I made some modifications (multiplied by -1 and add a fake phase plate for CTF) for that however the result is still strange.
Only because they are missing CTF info? No CTF is applied to the simulation so it doesn’t make sense to add it
this was written to test some new features :)
Author
Strange in what way?
The result of ReconstructParticleTomo doesn't resemble a ribosome structure.
Author
Ah, you will need to invert the rotations then
Sent from mobile - apologies for brevity
… On 21 Aug 2022, at 13:42, Xhark ***@***.***> wrote:
***@***.*** commented on this gist.
The result of ReconstructParticleTomo doesn't resemble a ribosome structure.
—
Reply to this email directly, view it on GitHub, or unsubscribe.
You are receiving this because you authored the thread.
Still something wrong.
Ah, you will need to invert the rotations then Sent from mobile - apologies for brevity
…
On 21 Aug 2022, at 13:42, Xhark @.> wrote: @. commented on this gist. The result of ReconstructParticleTomo doesn't resemble a ribosome structure. — Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you authored the thread.
Author
had a look - model_orientations is multiplying coordinates which are zyx not xyz, that's the source of the error
Does that mean rlnTomoProjX should be row 3 and rlnTomoProjZ should be row 1 of transformations?
had a look -
model_orientationsis multiplying coordinates which are zyx not xyz, that's the source of the error
Author
No, rlnTomoProjXYZW are fine, you would need to exchange columns 1 and 3 in the rotation matrices for the particles then exchange rows 1 and 3
something like below? looks like the transpose of the rotation matrix inside the affine transformation matrix.
inverse_transformations = transformations.copy()
for no in range(transformations.shape[0]):
inverse_transformations[no, :3, [2, 0]] = transformations[no, :3, [0, 2]]
inverse_transformations[no, [2, 0], :3] = inverse_transformations[no, [0, 2], :3]
No, rlnTomoProjXYZW are fine, you would need to exchange columns 1 and 3 in the rotation matrices for the particles then exchange rows 1 and 3
Author
I don't think so, here is some guidance...
transformationsare the (4, 4) projection matrices which related 3D position to 2D position in tilt-imagemodel_orientationsis a stack of (3, 3) rotation matrices which rotate the model (ribosome Ca coordinates) around its center
The problem is that I didn't need to keep track of the orientations carefully, I just wrote them out into the file for compatibility. To be sure they are correct you need to track what these matrices are transforming (xyz column vectors? zyx column vectors?) and how that data is placed into 3D/2D (axes of images are ordered zyx, slowest to fastest in memory)
All of these must be consistent and the conversion from rotation matrix to euler angles must follow the RELION convention for your experiment to produce the result you expect :) be careful, active rotation of coordinates is the inverse of the rotation of the coordinates of the fourier slice which is what relion does internally with its euler angles
Thanks for the explanation, hence it should be something below,
inverse_model_orientations = model_orientations.copy()
inverse_model_orientations[:, :, :, [2, 0]] = model_orientations[:, :, :, [0, 2]]
inverse_model_orientations[:, :, [2, 0], :] = inverse_model_orientations[:, :, [0, 2], :]
I don't think so, here is some guidance...
transformationsare the (4, 4) projection matrices which related 3D position to 2D position in tilt-imagemodel_orientationsis a stack of (3, 3) rotation matrices which rotate the model (ribosome Ca coordinates) around its centerThe problem is that I didn't need to keep track of the orientations carefully, I just wrote them out into the file for compatibility. To be sure they are correct you need to track what these matrices are transforming (xyz column vectors? zyx column vectors?) and how that data is placed into 3D/2D (axes of images are ordered zyx, slowest to fastest in memory)
All of these must be consistent and the conversion from rotation matrix to euler angles must follow the RELION convention for your experiment to produce the result you expect :) be careful, active rotation of coordinates is the inverse of the rotation of the coordinates of the fourier slice which is what relion does internally with its euler angles
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
The generated files can not be utilized by relion 4.1