{ "cells": [ { "cell_type": "markdown", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "How to calculate Global Elasticity?\n", "===================================\n", "\n", "To calculate global elasticity, ``dnaMD globalElasticity`` command can be used. However, it takes HDF5 file as input.\n", "Following steps can be performed to generate HDF5 files. The tutorial file can be downloaded here. We will prepare HDF5\n", "file for both free and bound DNA.\n", "\n", "Calculate stretching twisting and bending motions\n", "-------------------------------------------------\n", "\n", "Both free and bound DNA is superimposed on to the same DNA structure. Careful that bending calculation is fitting dependent.\n", "Therefore, at first we aligned both free and bound DNA to a common DNA structure as follows" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%bash\n", "\n", "# Align bound DNA\n", "echo 6 0 | gmx trjconv -f inputs/F1_complex_DNA.xtc -s inputs/dna.tpr -n inputs/dna.ndx -o complex_dna_aligned.xtc -fit rot+trans\n", "\n", "# Align free DNA\n", "echo 6 0 | gmx trjconv -f inputs/F1_free_DNA.xtc -s inputs/dna.tpr -n inputs/dna.ndx -o free_dna_aligned.xtc -fit rot+trans" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**1.** Run do_x3dna on DNA trajectory, ``-nofit`` is used because DNA is already superimposed to a common DNA structure using trjconv." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "raw_mimetype": "text/restructuredtext" }, "outputs": [], "source": [ "%%bash\n", "\n", "# For free DNA\n", "echo 2 | do_x3dna -f free_dna_aligned.xtc -s inputs/dna.tpr -n inputs/dna.ndx -ref -noavg -nofit -name free\n", "mv *_free.dat outputs/.\n", "\n", "# For bound DNA\n", "echo 2 | do_x3dna -f complex_dna_aligned.xtc -s inputs/dna.tpr -n inputs/dna.ndx -ref -noavg -nofit -name bound\n", "mv *_bound.dat outputs/." ] }, { "cell_type": "markdown", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "**2.** Run dnaMD to extract the parameters from do_x3dna output files and save as HDF5 file. Also calculate global axis,\n", "curvature and tangents." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "Loading parameters: ['shift', 'slide', 'rise', 'tilt', 'roll', 'twist']\n", "Reading file : L-BPS_free.dat\n", "Reading frame 5000\n", "Finished reading.... Total number of frame read = 5001\n", "\n", "Loading parameters: ['x-disp', 'y-disp', 'h-rise', 'inclination', 'tip', 'h-twist']\n", "Reading file : L-BPH_free.dat\n", "Reading frame 5000\n", "Finished reading.... Total number of frame read = 5001\n", "\n", "Loading parameters: ['helical x-axis', 'helical y-axis', 'helical z-axis']\n", "Reading file : HelAxis_free.dat\n", "Reading frame 5000\n", "Finished reading.... Total number of frame read = 5001\n", "Fitting spline curve on helical axis of frame 5000 out of 5001 frames\n", "Finished spline curve fitting...\n", " ... Calculating curvature and tangents \n", " ... Finished\n", "\n", "Loading parameters: ['shift', 'slide', 'rise', 'tilt', 'roll', 'twist']\n", "Reading file : L-BPS_bound.dat\n", "Reading frame 10000\n", "Finished reading.... Total number of frame read = 12001\n", "\n", "Loading parameters: ['x-disp', 'y-disp', 'h-rise', 'inclination', 'tip', 'h-twist']\n", "Reading file : L-BPH_bound.dat\n", "Reading frame 10000\n", "Finished reading.... Total number of frame read = 12001\n", "\n", "Loading parameters: ['helical x-axis', 'helical y-axis', 'helical z-axis']\n", "Reading file : HelAxis_bound.dat\n", "Reading frame 10000\n", "Finished reading.... Total number of frame read = 12001\n", "Fitting spline curve on helical axis of frame 12000 out of 12001 frames\n", "Finished spline curve fitting...\n", " ... Calculating curvature and tangents \n", " ... Finished\n" ] } ], "source": [ "%%bash\n", "\n", "# For free DNA\n", "dnaMD saveAsH5 -tbp 27 -i outputs/L-BPS_free.dat,outputs/L-BPH_free.dat,outputs/HelAxis_free.dat -o free_dna.h5\n", "dnaMD axisCurv -tbp 27 -bs 2 -be 25 -ctan -scp 100 -s 1000 -cta 30 -io free_dna.h5 -ap free_dna_axis.pdb\n", "\n", "# For bound DNA\n", "dnaMD saveAsH5 -tbp 27 -i outputs/L-BPS_bound.dat,outputs/L-BPH_bound.dat,outputs/HelAxis_bound.dat -o bound_dna.h5\n", "dnaMD axisCurv -tbp 27 -bs 2 -be 25 -ctan -scp 100 -s 1000 -cta 30 -io bound_dna.h5 -ap bound_dna_axis.pdb" ] }, { "cell_type": "markdown", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "Now, we have HDF5 files of both free and bounds DNA. It can be used for the calculation of elastic properties. These files\n", "can be used with either [dnaMD Python module](http://do-x3dna.readthedocs.io/en/latest/notebooks/calculate_elasticity_tutorial.html) or ``dnaMD globalElasticity``.\n", "\n", "\n", "Bending Stretching Twist Modulus\n", "--------------------------------\n", "\n", "Following command calculate Bending Stretching Twist modulus matrix. Output matrix will be stored in csv file.\n", "Elastic modulus matrix is printed as output and average values of contour length and cumulative twist angle\n", "is also printed." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "=========== Elastic Modulus Matrix ===========\n", "\n", " 337.233 28.506 -6.984 -20.710\n", " 28.506 373.863 7.977 39.509\n", " -6.984 7.977 1080.722 -96.198\n", " -20.710 39.509 -96.198 448.349\n", "\n", "=========== ====================== ===========\n", "\n", "========================= Average Values ==========================\n", "\n", "Bending-1 Angle Bending-2 Angle Contour Length Sum. Twist \n", " 0.021 0.068 5.358 9.543 \n", "\n", "=========== ====================== ====================== ===========\n" ] } ], "source": [ "%%bash\n", "\n", "dnaMD globalElasticity -i free_dna.h5 -tbp 27 -bs 4 -be 20 -estype BST -paxis X -o elastic_modulus_BST.csv" ] }, { "cell_type": "markdown", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "The above modulus matrix is in this form:\n", "\n", "$$\\text{modulus matrix} =\n", "\\begin{bmatrix}\n", "M_{Bx} & M_{Bx,By} & M_{Bx,S} & M_{Bx,T} \\\\\n", "M_{Bx,By} & M_{By} & M_{By,S} & M_{By,T} \\\\\n", "M_{Bx,S} & M_{By,S} & M_{S} & M_{S,T} \\\\\n", "M_{Bx,T} & M_{Bx,T} & M_{S,T} & M_{T}\n", "\\end{bmatrix}\n", "$$\n", "\n", "Where:\n", "\n", "* $M_{Bx}$ - Bending-1 stiffness in one plane\n", "* $M_{By}$ - Bending-2 stiffness in another orthogonal plane\n", "* $M_{S}$ - Stretch Modulus\n", "* $M_{T}$ - Twist rigidity\n", "* $M_{Bx,By}$ - Bending-1 and Bending-2 coupling\n", "* $M_{By,S}$ - Bending-2 and stretching coupling\n", "* $M_{S,T}$ - Stretching Twsiting coupling\n", "* $M_{Bx,S}$ - Bending-1 Stretching coupling\n", "* $M_{By,T}$ - Bending-2 Twisting coupling\n", "* $M_{Bx,T}$ - Bending-1 Twisting coupling\n", "\n", "\n", "Stretching Twist Modulus\n", "------------------------\n", "\n", "Following command calculate Bending Stretching Twist modulus matrix. Output matrix will be stored in csv file.\n", "Elastic modulus matrix is printed as output and average values of contour length and cumulative twist angle\n", "is also printed." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "=========== Elastic Modulus Matrix ===========\n", "\n", "1080.380 -97.578\n", "-97.578 442.494\n", "\n", "=========== ====================== ===========\n", "\n", "========================= Average Values ==========================\n", "\n", "Contour Length Sum. Twist \n", " 5.358 9.543 \n", "\n", "=========== ====================== ====================== ===========\n" ] } ], "source": [ "%%bash\n", "\n", "dnaMD globalElasticity -i free_dna.h5 -tbp 27 -bs 4 -be 20 -estype ST -o elastic_modulus_ST.csv" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The above modulus matrix is in this form:\n", "\n", "$$\\text{modulus matrix} =\n", "\\begin{bmatrix}\n", "M_{S} & M_{S,T} \\\\\n", "M_{S,T} & M_{T}\n", "\\end{bmatrix}\n", "$$\n", "\n", "Where:\n", "\n", "* $M_{S}$ - Stretch Modulus\n", "* $M_{T}$ - Twist rigidity\n", "* $M_{S,T}$ - Stretching Twsiting coupling\n", "\n", "Convergence in Modulus\n", "----------------------\n", "\n", "Same command can be used to calculate elasticity as a function of time with option \n", "``-ot/--output-time`` and save it in csv format file. This result can beused to check\n", "their convergence. \n", "\n", "* If this option is used, ``-fgap/--frame-gap`` is **an essential option**.\n", "* This options also gives final value and error as the output on display.\n", "* The output file is in **csv** format and can be opened as spreadsheet.\n", "\n", "**NOTE:** Elastic properties cannot be calculated using a single frame because \n", "fluctuations are required. Therefore, here time means trajectory between zero \n", "time to given time." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "==============================================\n", "Elasticity Value Error\n", "----------------------------------------------\n", "bend-1 337.233 3.932\n", "bend-2 373.863 4.935\n", "stretch 1080.722 24.015\n", "twist 448.349 29.356\n", "bend-1-bend-2 28.506 4.685\n", "bend-2-stretch 7.977 7.251\n", "stretch-twist -96.198 13.450\n", "bend-1-stretch -6.984 7.131\n", "bend-2-twist 39.509 35.576\n", "bend-1-twist -20.710 5.081\n", "==============================================\n", "\n", "=====================================\n", "Elastic modulus as a function of time\n", "=====================================\n", "#Time, bend-1, bend-2, stretch, twist, bend-1-bend-2, bend-2-stretch, stretch-twist, bend-1-stretch, bend-2-twist, bend-1-twist\n", "1000.000, 584.19567, 363.51886, 1121.88830, 1506.57584, 41.49043, -101.57973, 112.06246, -254.22188, -74.92882, 252.68516\n", "2000.000, 354.70978, 367.93829, 1032.49032, 662.58414, -18.43006, -53.55808, 9.82719, -81.25370, -207.81222, 111.70108\n", "3000.000, 380.26015, 359.63694, 967.63012, 543.68667, 5.70038, -76.61296, -8.92719, -88.41420, -95.81366, 56.24528\n", ".\n", ".\n", ".\n", "118000.000, 335.33614, 371.13759, 1082.41850, 449.84105, 27.64596, 10.92246, -97.17016, -7.92023, 36.92839, -17.34603\n", "119000.000, 335.79993, 373.00012, 1081.59948, 449.85747, 28.30290, 9.58966, -96.91853, -6.18165, 37.45497, -19.36140\n", "120000.000, 337.23300, 373.86276, 1080.72201, 448.34909, 28.50566, 7.97658, -96.19760, -6.98413, 39.50860, -20.71008\n" ] } ], "source": [ "%%bash\n", "\n", "dnaMD globalElasticity -i free_dna.h5 -tbp 27 -bs 4 -be 20 -estype BST -paxis X -fgap 100 -em block -ot modulus_time_BST.csv\n", "\n", "# Print first and last 3 line of output file\n", "echo \"=====================================\"\n", "echo \"Elastic modulus as a function of time\"\n", "echo \"=====================================\"\n", "head -4 modulus_time_BST.csv\n", "printf \".\\n.\\n.\\n\"\n", "tail -3 modulus_time_BST.csv" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Some of the plots from above data can be **found** [here](http://do-x3dna.readthedocs.io/en/latest/notebooks/calculate_elasticity_tutorial.html#Convergence-in-bending,-stretching-and-twisting-with-their-couplings)\n", "\n", "***\n", "\n", "Same as above but only for stretching and twisting motions." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "==============================================\n", "Elasticity Value Error\n", "----------------------------------------------\n", "stretch 1080.380 22.847\n", "twist 442.494 19.169\n", "stretch-twist -97.578 14.389\n", "==============================================\n", "\n", "=====================================\n", "Elastic modulus as a function of time\n", "=====================================\n", "#Time, stretch, twist, stretch-twist\n", "1000.000, 991.91168, 1373.35825, 200.50841\n", "2000.000, 1004.78021, 516.21083, 3.60912\n", "3000.000, 931.30826, 509.38341, -16.08887\n", ".\n", ".\n", ".\n", "118000.000, 1081.86818, 444.95480, -98.78361\n", "119000.000, 1081.20986, 444.61885, -98.34127\n", "120000.000, 1080.37979, 442.49438, -97.57808\n" ] } ], "source": [ "%%bash\n", "\n", "dnaMD globalElasticity -i free_dna.h5 -tbp 27 -bs 4 -be 20 -estype ST -paxis X -fgap 100 -gt \"gmx analyze\" -em \"block\" -ot modulus_time_ST.csv\n", "\n", "# Print first and last 3 line of output file\n", "echo \"=====================================\"\n", "echo \"Elastic modulus as a function of time\"\n", "echo \"=====================================\"\n", "head -4 modulus_time_ST.csv\n", "printf \".\\n.\\n.\\n\"\n", "tail -3 modulus_time_ST.csv" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Global deformation free energy\n", "------------------------------\n", "\n", "To caluclate global deformation enrgy, ``dnaMD globalEnergy`` can be used. At first, elastic matrix from reference\n", "DNA (most often free or unbound DNA) is calculated and subsequently this matrix is used to calculate deformation free\n", "energy of probe DNA (most often bound DNA).\n", "\n", "The deformation free energy is calculated using elastic matrix as follows\n", "\n", "$$G = \\frac{1}{2L_0}\\mathbf{xKx^T}$$\n", "\n", "$$\\mathbf{x} = \\begin{bmatrix}\n", " (\\theta^{x} - \\theta^{x}_0) & (\\theta^{y} - \\theta^{y}_0) & (L - L_0) & (\\phi - \\phi_0)\n", " \\end{bmatrix}$$\n", "\n", "Where, $\\mathbf{K}$, $\\theta^{x}_0$, $\\theta^{y}_0$, $L_0$ and $\\phi_0$ is calculated from reference DNA while $\\theta^{x}$, $\\theta^{y}$, $L$ and $\\phi$ is calculated for probe DNA from each frame.\n", "\n", "This command gives output energy as a function of time in csv file and also average energies with error." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "==============================================\n", "Energy(kJ/mol) Average Error\n", "----------------------------------------------\n", "full 27.714 0.724\n", "diag 27.012 0.751\n", "stretch 2.468 0.157\n", "twist 15.061 0.478\n", "st_coupling -0.743 0.040\n", "b1 1.477 0.180\n", "b2 8.005 0.266\n", "bend 9.718 0.386\n", "bs_coupling 0.034 0.002\n", "bt_coupling 0.825 0.020\n", "bb_coupling 0.235 0.020\n", "st 26.269 0.713\n", "bs 11.984 0.499\n", "bt 25.369 0.634\n", "==============================================\n", "\n", "===========================================================\n", "Deformation free energy of bound DNA as a function of time\n", "===========================================================\n", "#Time, full, diag, stretch, twist, st_coupling, b1, b2, bend, bs_coupling, bt_coupling, bb_coupling, st, bs, bt\n", "0.000, 29.76525, 29.12752, 3.05586, 16.01787, -0.96687, 0.27243, 9.78136, 10.18484, 0.05805, 1.09664, 0.13105, 28.16064, 13.16770, 27.16829\n", "10.000, 27.69000, 26.82145, 1.78113, 13.65491, -0.68154, 3.00373, 8.38168, 11.78823, 0.02173, 0.69127, 0.40282, 26.13991, 13.18827, 25.73159\n", "20.000, 32.65326, 31.73848, 3.36328, 15.22752, -0.98900, 0.93046, 12.21722, 13.41835, 0.05997, 1.11574, 0.27067, 30.74949, 16.57094, 29.49094\n", ".\n", ".\n", ".\n", "49980.000, 28.48093, 28.59762, 4.40662, 16.38502, -1.17429, 0.45068, 7.35530, 7.95214, 0.05514, 0.91465, 0.14616, 27.42332, 12.26774, 25.10564\n", "49990.000, 33.18350, 32.91519, 3.56658, 21.03477, -1.19700, 0.05849, 8.25534, 8.36962, 0.06281, 1.21256, 0.05579, 31.71819, 11.94323, 30.56117\n", "50000.000, 30.81292, 29.94787, 1.86541, 19.01796, -0.82313, 0.67654, 8.38796, 9.25575, 0.03664, 1.02777, 0.19124, 29.12474, 10.96655, 29.11023\n" ] } ], "source": [ "%%bash\n", "\n", "dnaMD globalEnergy -ir free_dna.h5 -ip bound_dna.h5 -tbp 27 -bs 4 -be 20 -estype BST -et \"all\" -paxis X -gt \"gmx analyze\" -em \"block\" -o energy_all_BST.csv\n", "\n", "# Print first and last 3 line of output file\n", "echo \"===========================================================\"\n", "echo \"Deformation free energy of bound DNA as a function of time\"\n", "echo \"===========================================================\"\n", "head -4 energy_all_BST.csv\n", "printf \".\\n.\\n.\\n\"\n", "tail -3 energy_all_BST.csv" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Some of the plots from above data can be **found** [here](http://do-x3dna.readthedocs.io/en/latest/notebooks/calculate_elasticity_tutorial.html#Deformation-free-energy-of-bound-DNA)\n", "\n", "***\n", "\n", "Deformation free energy can be calculated as the following terms:\n", "\n", "* ``full`` : Use entire elastic matrix -- all motions with their coupling\n", "* ``diag`` : Use diagonal of elastic matrix -- all motions but no coupling\n", "* ``b1`` : Only bending-1 motion\n", "* ``b2`` : Only bending-2 motion\n", "* ``stretch`` : Only stretching motion\n", "* ``twist`` : Only Twisting motions\n", "* ``st_coupling`` : Only stretch-twist coupling motion\n", "* ``bs_coupling`` : Only Bending and stretching coupling\n", "* ``bt_coupling`` : Only Bending and Twisting coupling\n", "* ``bb_coupling`` : Only bending-1 and bending-2 coupling\n", "* ``bend`` : Both bending motions with their coupling\n", "* ``st`` : Stretching and twisting motions with their coupling\n", "* ``bs`` : Bending (b1, b2) and stretching motions with their coupling\n", "* ``bt`` : Bending (b1, b2) and twisting motions with their coupling\n", "\n", "When ``all`` is used, all above terms were calculated." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "==============================================\n", "Energy(kJ/mol) Average Error\n", "----------------------------------------------\n", "full 15.398 0.353\n", "diag 17.211 0.408\n", "stretch 2.581 0.137\n", "twist 14.630 0.362\n", "st_coupling -0.906 0.037\n", "==============================================\n", "\n", "===========================================================\n", "Deformation free energy of bound DNA as a function of time\n", "===========================================================\n", "#Time, full, diag, stretch, twist, st_coupling\n", "0.000, 24.65120, 27.93891, 4.01772, 23.92119, -1.64386\n", "10.000, 16.45031, 18.42131, 2.11880, 16.30251, -0.98550\n", "20.000, 22.50793, 25.46474, 3.54666, 21.91808, -1.47841\n", ".\n", ".\n", ".\n", "49980.000, 13.37563, 15.94228, 5.74292, 10.19936, -1.28332\n", "49990.000, 17.67070, 20.02852, 2.88301, 17.14551, -1.17891\n", "50000.000, 16.73218, 18.89026, 2.53138, 16.35888, -1.07904\n" ] } ], "source": [ "%%bash\n", "\n", "dnaMD globalEnergy -ir free_dna.h5 -ip bound_dna.h5 -tbp 27 -bs 4 -be 24 -estype ST -et \"all\" -gt \"gmx analyze\" -em \"block\" -o energy_all_ST.csv\n", "\n", "# Print first and last 3 line of output file\n", "echo \"===========================================================\"\n", "echo \"Deformation free energy of bound DNA as a function of time\"\n", "echo \"===========================================================\"\n", "head -4 energy_all_ST.csv\n", "printf \".\\n.\\n.\\n\"\n", "tail -3 energy_all_ST.csv" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Same as above but only for stretching and twisting motions.\n", "\n", "Local elastic properties\n", "------------------------\n", "\n", "Local elastic properties can be caluclated using either local base-step parameters or local helical base-step parameters.\n", "\n", "In case of **base-step** parameters: Shift ($Dx$), Slide ($Dy$), Rise ($Dz$), Tilt ($\\tau$), Roll ($\\rho$) and Twist ($\\omega$), following elastic matrix is calculated.\n", "\n", "$$\n", "\\mathbf{K}_{base-step} = \\begin{bmatrix}\n", "K_{Dx} & K_{Dx,Dy} & K_{Dx,Dz} & K_{Dx,\\tau} & K_{Dx,\\rho} & K_{Dx,\\omega} \\\\\n", "K_{Dx,Dy} & K_{Dy} & K_{Dy,Dz} & K_{Dy,\\tau} & K_{Dy,\\rho} & K_{Dy,\\omega} \\\\\n", "K_{Dx,Dz} & K_{Dy,Dz} & K_{Dz} & K_{Dz,\\tau} & K_{Dz,\\rho} & K_{Dz,\\omega} \\\\\n", "K_{Dx,\\tau} & K_{Dy,\\tau} & K_{Dz,\\tau} & K_{\\tau} & K_{\\tau, \\rho} & K_{\\tau,\\omega} \\\\\n", "K_{Dx,\\rho} & K_{Dy,\\rho} & K_{Dz,\\rho} & K_{\\tau, \\rho} & K_{\\rho} & K_{\\rho,\\omega} \\\\\n", "K_{Dx,\\omega} & K_{Dy,\\omega} & K_{Dz,\\omega} & K_{\\tau, \\omega} & K_{\\rho, \\omega} & K_{\\omega} \\\\\n", "\\end{bmatrix}\n", "$$\n", "\n", "In case of **helical-base-step** parameters: x-displacement ($dx$), y-displacement ($dy$), h-rise ($h$), inclination ($\\eta$), tip ($\\theta$) and twist ($\\Omega$), following elastic matrix is calculated.\n", "\n", "$$\n", "\\mathbf{K}_{helical-base-step} = \\begin{bmatrix}\n", "K_{dx} & K_{dx,dy} & K_{dx,h} & K_{dx,\\eta} & K_{dx,\\theta} & K_{dx,\\Omega} \\\\\n", "K_{dx,dy} & K_{dy} & K_{dy,h} & K_{dy,\\eta} & K_{dy,\\theta} & K_{dy,\\Omega} \\\\\n", "K_{dx,h} & K_{dy,h} & K_{h} & K_{h,\\eta} & K_{h,\\theta} & K_{h,\\Omega} \\\\\n", "K_{dx,\\eta} & K_{dy,\\eta} & K_{h,\\eta} & K_{\\eta} & K_{\\eta, \\theta} & K_{\\eta,\\Omega} \\\\\n", "K_{dx,\\theta} & K_{dy,\\theta} & K_{h,\\theta} & K_{\\eta, \\theta} & K_{\\theta} & K_{\\theta,\\Omega} \\\\\n", "K_{dx,\\Omega} & K_{dy,\\Omega} & K_{h,\\Omega} & K_{\\eta, \\Omega} & K_{\\theta, \\Omega} & K_{\\Omega} \\\\\n", "\\end{bmatrix}\n", "$$" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "=========== ============== Elastic Matrix =============== ===========\n", "\n", " 150.32688 -12.74081 -4.99484 -0.86502 -0.08939 -0.16720\n", " -12.74081 129.99578 30.42056 0.27454 -0.58620 -2.06912\n", " -4.99484 30.42056 504.89222 0.65083 -0.06531 0.21751\n", " -0.86502 0.27454 0.65083 0.03433 -0.00104 0.00052\n", " -0.08939 -0.58620 -0.06531 -0.00104 0.01295 0.01427\n", " -0.16720 -2.06912 0.21751 0.00052 0.01427 0.05907\n", "\n", "=========== ====================== ====================== ===========\n", "\n" ] } ], "source": [ "%%bash\n", "\n", "dnaMD localElasticity -i free_dna.h5 -tbp 27 -bs 4 -be 7 -o local_elasticity_4-7bps.csv" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Local elastic properties as a function of time\n", "----------------------------------------------\n", "\n", "Same command can be used to calculate elasticity as a function of time with option \n", "``-ot/--output-time`` and save it in csv format file. This result can be used to check\n", "their convergence. \n", "\n", "* If this option is used, ``-fgap/--frame-gap`` is **an essential option**.\n", "* The output file is in **csv** format and can be opened as spreadsheet.\n", "\n", "**NOTE:** Elastic properties cannot be calculated using a single frame because \n", "fluctuations are required. Therefore, here time means trajectory between zero \n", "time to given time." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "=====================================\n", "Elastic matrix as a function of time\n", "=====================================\n", "#Time, shift, slide, rise, tilt, roll, twist, ... ... ...\n", "2000.000, 158.96347, 127.71194, 503.26390, 0.03143, 0.01542, 0.05937, ... ... ...\n", "4000.000, 150.93805, 126.16929, 467.98254, 0.03050, 0.01409, 0.06100, ... ... ...\n", "6000.000, 147.43936, 132.39932, 480.47463, 0.03256, 0.01358, 0.06143, ... ... ...\n", ".\n", ".\n", ".\n", "116000.000, 151.11253, 129.64049, 504.93262, 0.03417, 0.01293, 0.05897, ... ... ...\n", "118000.000, 150.77496, 129.62910, 504.85789, 0.03426, 0.01296, 0.05910, ... ... ...\n", "120000.000, 150.32688, 129.99578, 504.89222, 0.03433, 0.01295, 0.05907, ... ... ...\n" ] } ], "source": [ "%%bash\n", "\n", "dnaMD localElasticity -i free_dna.h5 -tbp 27 -bs 4 -be 7 -fgap 200 -ot local_elasticity_time_4-7bps.csv\n", "\n", "# Print first 3 and last 3 line of first seven columns from output file\n", "echo \"=====================================\"\n", "echo \"Elastic matrix as a function of time\"\n", "echo \"=====================================\"\n", "head -4 local_elasticity_time_4-8bps.csv | awk '{print $1, $2, $3, $4, $5, $6 ,$7, \"... ... ...\"}'\n", "printf \".\\n.\\n.\\n\"\n", "tail -3 local_elasticity_time_4-8bps.csv | awk '{print $1, $2, $3, $4, $5, $6 ,$7, \"... ... ...\"}'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Local elastic properties of the consecutive overlapped DNA segments\n", "--------------------------------------------------------------------\n", "\n", "Above method gives local elasticities of a small local segment of the DNA. However, we mostly interested\n", "in large segment of the DNA. This large segment can be further divided into smaller local segments. \n", "For these smaller segments local elasticities can be calculated. Here these segments overlapped with each other.\n", "\n", "Same command can be used to calculate elasticity as a function of time with option \n", "``-os/--output-segments`` and save it in csv format file.\n", "\n", "* If this option is used, ``-fgap/--frame-gap`` is **an essential option**.\n", "* The output file is in **csv** format and can be opened as spreadsheet.\n" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "========================================\n", "Elasticity as a function of DNA segments\n", "========================================\n", "#bps, shift, shift-error, slide, slide-error, rise, rise-error, ... ... ...\n", "4-7, 150.32688, 0.25654, 129.99578, 0.27262, 504.89222, 1.61893, ... ... ...\n", "5-8, 162.67911, 0.46671, 122.65855, 0.47513, 510.18741, 1.36958, ... ... ...\n", "6-9, 172.50587, 2.36981, 119.84961, 1.46698, 510.68380, 2.57983, ... ... ...\n", ".\n", ".\n", ".\n", "15-18, 153.98831, 2.02295, 137.21947, 0.29044, 546.15664, 1.31991, ... ... ...\n", "16-19, 188.51349, 1.91188, 148.97261, 0.46969, 521.74643, 1.14698, ... ... ...\n", "17-20, 174.28796, 2.07481, 134.04039, 0.67123, 573.20011, 2.54414, ... ... ...\n" ] } ], "source": [ "%%bash\n", "\n", "dnaMD localElasticity -i free_dna.h5 -tbp 27 -bs 4 -be 20 -span 4 -fgap 200 -gt \"gmx analyze\" -em \"acf\" -os local_elasticity_segments.csv\n", "\n", "# Print first 3 and last 3 line of first seven columns from output file\n", "echo \"========================================\"\n", "echo \"Elasticity as a function of DNA segments\"\n", "echo \"========================================\"\n", "head -4 local_elasticity_segments.csv | awk '{print $1, $2, $3, $4, $5, $6 ,$7, \"... ... ...\"}'\n", "printf \".\\n.\\n.\\n\"\n", "tail -3 local_elasticity_segments.csv | awk '{print $1, $2, $3, $4, $5, $6 ,$7, \"... ... ...\"}'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Local deformation energy of a local small segment\n", "-------------------------------------------------\n", "At first, elastic matrix from reference DNA (most often free or unbound DNA) is calculated \n", "and subsequently this matrix is used to calculate deformation free energy of probe DNA \n", "(most often bound DNA).\n", "\n", "$$G = \\frac{1}{2}\\mathbf{xKx^T}$$\n", "\n", "When ``helical='False'``\n", "\n", "$$\\mathbf{K} = \\mathbf{K}_{base-step}$$\n", "\n", "$$\\mathbf{x} = \\begin{bmatrix}\n", " (Dx_{i}-Dx_0) & (Dy_i - Dy_0) & (Dz_i - Dz_0) & (\\tau_i - \\tau_0) &\n", " (\\rho_i - \\rho_0) & (\\omega_i - \\omega_0)\n", " \\end{bmatrix}$$\n", "\n", "\n", "When ``helical='True'``\n", "\n", "$$\\mathbf{K} = \\mathbf{K}_{helical-base-step}$$\n", "\n", "$$\\mathbf{x} = \\begin{bmatrix}\n", " (dx_{i}-dx_0) & (dy_i - dy_0) & (h_i - h_0) & (\\eta_i - \\eta_0) &\n", " (\\theta_i - \\theta_0) & (\\Omega_i - \\Omega_0)\n", " \\end{bmatrix}$$\n" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "==============================================\n", "Energy(kJ/mol) Average Error\n", "----------------------------------------------\n", "full 19.143 0.671\n", "diag 43.865 1.839\n", "shift 1.366 0.064\n", "slide 15.994 1.069\n", "rise 1.194 0.030\n", "tilt 1.705 0.040\n", "roll 6.365 0.365\n", "twist 17.242 0.781\n", "==============================================\n", "\n", "===============================================\n", "Local deformation energy as a function of time\n", "===============================================\n", "#Time, full, diag, shift, slide, rise, tilt, roll, twist\n", "0.000, 32.99568, 97.62372, 0.85609, 43.00855, 0.55182, 0.05174, 8.64099, 44.51453\n", "10.000, 24.91292, 63.21362, 0.18680, 31.85144, 0.31550, 2.74114, 5.80446, 22.31428\n", "20.000, 34.19821, 78.73263, 0.77440, 31.54860, 4.31769, 0.19352, 8.74777, 33.15064\n", ".\n", ".\n", ".\n", "49980.000, 17.43224, 50.92180, 0.07917, 12.95852, 2.46789, 0.92035, 5.74851, 28.74736\n", "49990.000, 18.86625, 58.72850, 3.25713, 25.11722, 0.01483, 1.27080, 1.81182, 27.25671\n", "50000.000, 19.01918, 49.56285, 2.16952, 13.54590, 1.86029, 4.52988, 4.38336, 23.07391\n" ] } ], "source": [ "%%bash\n", "\n", "dnaMD localEnergy -ir free_dna.h5 -ip bound_dna.h5 -tbp 27 -bs 10 -be 13 -et all -gt \"gmx analyze\" -em \"block\" -o local_energy_time_4-7bps.csv\n", "\n", "# Print first 3 and last 3 line from output file\n", "echo \"===============================================\"\n", "echo \"Local deformation energy as a function of time\"\n", "echo \"===============================================\"\n", "head -4 local_energy_time_4-7bps.csv\n", "printf \".\\n.\\n.\\n\"\n", "tail -3 local_energy_time_4-7bps.csv" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Some of the plots from above data can be **found** \n", "[here](http://do-x3dna.readthedocs.io/en/latest/notebooks/calculate_elasticity_tutorial.html#Local-deformation-energy-of-a-local-small-segment)\n", "\n", "***\n", "\n", "Same as the above but energy is calculated using helical base-step parameters" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "==============================================\n", "Energy(kJ/mol) Average Error\n", "----------------------------------------------\n", "full 18.541 0.607\n", "diag 82.103 3.204\n", "x-disp 38.632 1.974\n", "y-disp 0.940 0.024\n", "h-rise 6.768 0.333\n", "inclination 15.242 0.857\n", "tip 1.288 0.030\n", "h-twist 19.232 0.825\n", "==============================================\n", "\n", "======================================================\n", "Local helical deformation energy as a function of time\n", "======================================================\n", "#Time, full, diag, shift, slide, rise, tilt, roll, twist\n", "0.000, 32.99568, 97.62372, 0.85609, 43.00855, 0.55182, 0.05174, 8.64099, 44.51453\n", "10.000, 24.91292, 63.21362, 0.18680, 31.85144, 0.31550, 2.74114, 5.80446, 22.31428\n", "20.000, 34.19821, 78.73263, 0.77440, 31.54860, 4.31769, 0.19352, 8.74777, 33.15064\n", ".\n", ".\n", ".\n", "49980.000, 17.43224, 50.92180, 0.07917, 12.95852, 2.46789, 0.92035, 5.74851, 28.74736\n", "49990.000, 18.86625, 58.72850, 3.25713, 25.11722, 0.01483, 1.27080, 1.81182, 27.25671\n", "50000.000, 19.01918, 49.56285, 2.16952, 13.54590, 1.86029, 4.52988, 4.38336, 23.07391\n" ] } ], "source": [ "%%bash\n", "\n", "dnaMD localEnergy -ir free_dna.h5 -ip bound_dna.h5 -tbp 27 -bs 10 -be 13 -et all -helical -gt \"gmx analyze\" -em \"block\" -o local_helical_energy_time_4-7bps.csv\n", "\n", "# Print first 3 and last 3 line from output file\n", "echo \"======================================================\"\n", "echo \"Local helical deformation energy as a function of time\"\n", "echo \"======================================================\"\n", "head -4 local_energy_time_4-7bps.csv\n", "printf \".\\n.\\n.\\n\"\n", "tail -3 local_energy_time_4-7bps.csv" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Deformation energy of the consecutive overlapped DNA segments\n", "-------------------------------------------------------------\n", "\n", "Above method gives energy of a small local segment of the DNA. \n", "However, we mostly interested in large segment of the DNA. This large segment \n", "can be further divided into smaller local segments. For these smaller segments \n", "local deformation energy can be calculated. Here these segments overlapped with each other." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "==================================================\n", "Local deformation energy as a function of segments\n", "==================================================\n", "#bps, full, full-error, diag, diag-error, shift, shift-error, slide, slide-error, rise, rise-error, tilt, tilt-error, roll, roll-error, twist, twist-error\n", "4-7, 10.44541, 0.15599, 16.11461, 0.50862, 0.95000, 0.02392, 8.97490, 0.44581, 1.17289, 0.13403, 1.71030, 0.05021, 1.08986, 0.03874, 2.21666, 0.17392\n", "5-8, 9.49626, 0.25053, 19.32846, 0.70809, 1.14060, 0.02737, 9.59726, 0.51038, 1.21808, 0.18759, 1.75363, 0.04521, 1.11454, 0.05955, 4.50435, 0.19880\n", "6-9, 10.52567, 0.28519, 20.09691, 0.54228, 1.22326, 0.02970, 11.45709, 0.43514, 1.51242, 0.08500, 1.90555, 0.04658, 1.31208, 0.11407, 2.68650, 0.09203\n", ".\n", ".\n", ".\n", "15-18, 13.98537, 0.52298, 41.05547, 2.63829, 1.04577, 0.02497, 17.01431, 0.86024, 2.44413, 0.14371, 1.78045, 0.04873, 1.78854, 0.20280, 16.98227, 2.09869\n", "16-19, 14.04847, 0.37844, 39.34525, 1.49491, 1.97148, 0.09511, 14.20378, 0.94662, 2.52573, 0.08194, 1.70940, 0.04099, 0.98835, 0.09072, 17.94652, 1.15043\n", "17-20, 7.89839, 0.32320, 14.91991, 0.78938, 1.03028, 0.02583, 5.81921, 0.49549, 1.40673, 0.06566, 1.46641, 0.03364, 1.05240, 0.14133, 4.14489, 0.38861\n" ] } ], "source": [ "%%bash\n", "\n", "dnaMD localEnergy -ir free_dna.h5 -ip bound_dna.h5 -tbp 27 -bs 4 -be 20 -span 4 -et all -gt \"gmx analyze\" -em \"block\" -os local_energy_segments.csv\n", "\n", "# Print first 3 and last 3 line from output file\n", "echo \"==================================================\"\n", "echo \"Local deformation energy as a function of segments\"\n", "echo \"==================================================\"\n", "head -4 local_energy_segments.csv\n", "printf \".\\n.\\n.\\n\"\n", "tail -3 local_energy_segments.csv" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Some of the plots from above data can be **found** \n", "[here](http://do-x3dna.readthedocs.io/en/latest/notebooks/calculate_elasticity_tutorial.html#Deformation-energy-of-the-consecutive-overlapped-DNA-segments)\n", "\n", "***\n", "\n", "\n", "Same as the above but energy is calculated using helical base-step parameters" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "==========================================================\n", "Local helical deformation energy as a function of segments\n", "==========================================================\n", "#bps, full, full-error, diag, diag-error, x-disp, x-disp-error, y-disp, y-disp-error, h-rise, h-rise-error, inclination, inclination-error, tip, tip-error, h-twist, h-twist-error\n", "4-7, 9.94039, 0.17649, 17.73868, 0.51652, 9.00892, 0.37199, 1.19946, 0.04468, 1.54174, 0.13701, 1.69527, 0.06446, 1.73424, 0.04823, 2.55906, 0.16107\n", "5-8, 9.01980, 0.22539, 19.49642, 0.57573, 8.19696, 0.26473, 1.36335, 0.06127, 1.48339, 0.27526, 1.60444, 0.08879, 1.76889, 0.04544, 5.07939, 0.32815\n", "6-9, 10.30985, 0.26888, 17.00383, 0.38032, 7.38157, 0.25950, 1.96202, 0.28106, 0.99682, 0.02458, 1.68845, 0.14166, 1.89859, 0.04606, 3.07639, 0.16278\n", ".\n", ".\n", ".\n", "15-18, 13.74199, 0.55739, 53.29447, 4.60218, 27.45388, 1.71961, 0.66820, 0.01620, 1.17571, 0.02924, 3.37921, 0.40760, 1.29946, 0.03306, 19.31802, 2.43431\n", "16-19, 12.60741, 0.33640, 46.34514, 2.01539, 19.89670, 0.90439, 0.90573, 0.04403, 1.01976, 0.02504, 1.68335, 0.09711, 1.12857, 0.06776, 21.71103, 1.22043\n", "17-20, 7.54442, 0.33087, 15.27288, 0.76005, 5.60088, 0.40791, 0.91852, 0.02377, 1.13807, 0.02707, 1.66290, 0.20744, 1.25180, 0.02884, 4.70071, 0.43934\n" ] } ], "source": [ "%%bash\n", "\n", "dnaMD localEnergy -ir free_dna.h5 -ip bound_dna.h5 -tbp 27 -bs 4 -be 20 -span 4 -et all -helical -gt \"gmx analyze\" -em \"block\" -os local_helical_energy_segments.csv\n", "\n", "# Print first 3 and last 3 line from output file\n", "echo \"==========================================================\"\n", "echo \"Local helical deformation energy as a function of segments\"\n", "echo \"==========================================================\"\n", "head -4 local_helical_energy_segments.csv\n", "printf \".\\n.\\n.\\n\"\n", "tail -3 local_helical_energy_segments.csv" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.11.11" } }, "nbformat": 4, "nbformat_minor": 4 }