{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modifying Solenoids\n",
    "\n",
    "Solenoids can be modified to change the density of the DNA.\n",
    "Here we explore how this works.\n",
    "\n",
    "The call to the Solenoid function is as follows:\n",
    "```\n",
    "dna.Solenoid(\n",
    "    voxelheight: float = 750,\n",
    "    radius: float = 100,\n",
    "    nhistones: int = 38,\n",
    "    histone_angle: float = 50,\n",
    "    twist: bool = False,\n",
    "    chain: int = 0,\n",
    ")\n",
    "```\n",
    "\n",
    "These numbers are a realistic model for how much Solenoidal DNA packs.\n",
    "\n",
    "Unchanged, this will build a solenoid\n",
    "- 750 nm tall\n",
    "- histones centered 100A from the central axis\n",
    "- 38 histones\n",
    "- histones tilted at 50°.\n",
    "\n",
    "Histones are always placed 60° apart going around the circle\n",
    "\n",
    "*It is recommended that you run this example inside a Jupyter environment*\n",
    "*rather than a VSCode or similar environment*\n",
    "\n",
    "This requires the mayavi jupyter extension\n",
    "`jupyter nbextension install --py mayavi --user`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "from pathlib import Path\n",
    "\n",
    "try:\n",
    "    from fractaldna.dna_models import dnachain as dna\n",
    "except (ImportError, ModuleNotFoundError):\n",
    "    sys.path.append(str(Path.cwd().parent.parent.parent))\n",
    "    from fractaldna.dna_models import dnachain as dna\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from mayavi import mlab\n",
    "\n",
    "# Disable this option for interactive rendering\n",
    "mlab.options.offscreen = True\n",
    "\n",
    "# Enable this option for an interactive notebook\n",
    "# mlab.init_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Varying Voxel Height\n",
    "Varying the voxel height will compact and extend the DNA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "chain_750 = dna.Solenoid(voxelheight=750)\n",
    "# MayaVI plots are best for visualisation here\n",
    "plot = chain_750.to_line_plot()\n",
    "\n",
    "# Save the figure\n",
    "plot.scene.save_jpg(\"single_solenoid_750A.jpg\")\n",
    "\n",
    "# In an interactive notebook, you can refer to the figure to\n",
    "# interact with it.\n",
    "plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "750 Å section, 38 histones\n",
    "\n",
    "![Single 750A solenoid](single_solenoid_750A.jpg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make a straight solenoid in a 500 Å box\n",
    "chain_500 = dna.Solenoid(voxelheight=500)\n",
    "# MayaVI plots are best for visualisation here\n",
    "plot = chain_500.to_line_plot()\n",
    "\n",
    "# Save the figure\n",
    "plot.scene.save_jpg(\"single_solenoid_500A.jpg\")\n",
    "\n",
    "# In an interactive notebook, you can refer to the figure to\n",
    "# interact with it.\n",
    "plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "500 Å section, 38 histones - the DNA is more compacted (and the DNA could have overlaps).\n",
    "\n",
    "![Shorter Solenoid](single_solenoid_500A.jpg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Changing the number of Histones\n",
    "When making DNA shorter (or longer), you should change the number of histones appropriately.\n",
    "\n",
    "In general, you'll want to keep the linear density of histones roughly constant (around 51/1000Å)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make a turned solenoid in a 750 Å box\n",
    "chain = dna.Solenoid(voxelheight=500, nhistones=25)\n",
    "# MayaVI plots are best for visualisation here\n",
    "plot = chain.to_line_plot()\n",
    "plot.scene.save_jpg(\"single_solenoid_500A_25histones.jpg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A 500Å long solenoid\n",
    "\n",
    "![500A Solenoid with 25 histones](single_solenoid_500A_25histones.jpg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*If you need to make DNA a little denser, try increasing the linear density of DNA*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make a turned solenoid in a 750 Å box\n",
    "chain = dna.Solenoid(voxelheight=750, nhistones=42)\n",
    "# MayaVI plots are best for visualisation here\n",
    "plot = chain.to_line_plot()\n",
    "plot.scene.save_jpg(\"single_solenoid_750A_42histones.jpg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "750Å strand, 42 histones\n",
    "\n",
    "![Straight Single Solenoid](single_solenoid_750A_42histones.jpg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Turned Solenoids follow the same principles\n",
    "\n",
    "They are actually based on the spacing that a straight solenoid would have,\n",
    "so the `voxelheight` and `nhistones` correspond to their straightened\n",
    "equivalents, and the linear density of histones is what is conserved.\n",
    "\n",
    "The radius of curvature is _half the voxel height_.\n",
    "\n",
    "Be particularly careful to ensure there are no overlapping molecules on the internal curve."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "chain = dna.TurnedSolenoid(voxelheight=750, nhistones=42)\n",
    "# MayaVI plots are best for visualisation here\n",
    "plot = chain.to_line_plot()\n",
    "plot.scene.save_jpg(\"turned_solenoid_750A_42histones.jpg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "750Å strand, 42 histones\n",
    "\n",
    "\n",
    "![Turned Solenoid 750A, 42 histones](turned_solenoid_750A_42histones.jpg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*It is important to inspect the result to make sure there aren't placements where you don't expect*\n",
    "\n",
    "With a voxelheight of 750Å, we expect bounds in the range\n",
    "\n",
    "-375 < X < 375\n",
    "-375 < Y < 375\n",
    "0 < Z < 750"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "chain_df = chain.to_frame()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "chain_df[\"max_extent\"] = chain_df[[\"size_x\", \"size_y\", \"size_z\"]].max(axis=1)\n",
    "\n",
    "# What are the bounds of the box?\n",
    "print(\n",
    "    \"X Bounds:\",\n",
    "    (chain_df[\"pos_x\"] - chain_df[\"max_extent\"]).min(),\n",
    "    (chain_df[\"pos_x\"] - chain_df[\"max_extent\"]).max(),\n",
    ")\n",
    "print(\n",
    "    \"Y Bounds:\",\n",
    "    (chain_df[\"pos_y\"] - chain_df[\"max_extent\"]).min(),\n",
    "    (chain_df[\"pos_y\"] - chain_df[\"max_extent\"]).max(),\n",
    ")\n",
    "print(\n",
    "    \"Z Bounds:\",\n",
    "    (chain_df[\"pos_z\"] - chain_df[\"max_extent\"]).min(),\n",
    "    (chain_df[\"pos_z\"] - chain_df[\"max_extent\"]).max(),\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Example Output\n",
    "```\n",
    "X Bounds: -130.5016837010315 371.3927605800866\n",
    "Y Bounds: -152.83468961083372 148.26998187713423\n",
    "Z Bounds: 0.04677445755690002 483.71674638609016\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exporting Solenoidal DNA to CSV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This too can be exported to a dataframe of basepairs\n",
    "\n",
    "# let's shift in so that it starts at z=-750/2, not at z=0\n",
    "chain.translate([0, 0, -750.0 / 2.0])\n",
    "\n",
    "# and then make it into a data frame\n",
    "chain.to_frame()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# And a second data frame of histones\n",
    "# This can be joined to the base pairs frame, assuming a sufficient handling of\n",
    "# the missing strand index, and the relation between the base pair and histone index.\n",
    "chain.histones_to_frame()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.13 ('fractaldna-ygFQyeT7-py3.9')",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  },
  "vscode": {
   "interpreter": {
    "hash": "4d48c20f365c2ac57f4f2737b0871661a9a1c46fb4f54286f41a4c7d06c7cec5"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}