Skip to content

Commit

Permalink
elasticity tools and tutorials added
Browse files Browse the repository at this point in the history
  • Loading branch information
@rjdkmr
rjdkmr committed Oct 17, 2018
1 parent 83cf312 commit b43e14d
Show file tree
Hide file tree
Showing 25 changed files with 3,871 additions and 665 deletions.
13 changes: 12 additions & 1 deletion README.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,12 +33,23 @@ It consists of three main components:
**users with programming experiences**.


**Last Update: Jan. 2018**
**Last Update: Oct. 2018**

**For detailed documentation about the do_x3dna, please visit:** |do_x3dna homepage|.

**For Questions and Discussions, please visit:** `do_x3dna forum <https://groups.google.com/forum/#!forum/do_x3dna>`_.

Release Note 2018
-----------------

* `Added new dnaMD tools <http://do-x3dna.readthedocs.io/en/latest/dnaMD_usage.html#commands-table>`_

* Added support to calculate both `global and local Elastic properties <http://do-x3dna.readthedocs.io/en/latest/about_dna_elasticity.html>`_

* Added tutorials to calculate elastic properties and deformation energy
using `dnaMD tool <http://do-x3dna.readthedocs.io/en/latest/global_elasticity.html>`_
and `dnaMD Python module <http://do-x3dna.readthedocs.io/en/latest/notebooks/calculate_elasticity_tutorial.html>`_.

Release Note 2017
-----------------

Expand Down
204 changes: 109 additions & 95 deletions dnaMD/dnaMD/commands/globalElasticity.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,28 +49,28 @@

description="""DESCRIPTION
===========
Parameter distribution during simulation
Global Elastic Properties of the DNA
This can be used to calculate distribution of a parameter of either a specific base-pair/step
or over a DNA segment during the simulations.
This can be used to Global Elastic Properties of the DNA from the simulations.
"""

inputFileHelp=\
"""Name of input file (from do_x3dna or hdf5 file).
This file should contain the required parameters. It can be a file either
produced from do_x3dna or hdf5 storage file.
"""Name of input file (hdf5 file).
This file should contain the required parameters. It should be hdf5 storage file.
"""

outputFileHelp=\
"""Name of output file.
The elasticity modulus will be written in this file.
"""Name of output file with csv extension.
This file will contain the elasticity modulus matrix where values will be separated
by comma. Since modulus matrix is also shown as screen output, this option is not
necessary.
"""

vsTimeHelp=\
"""Calculate moduli as a function of time.
"""Calculate elasticity as a function of time and save in this csv format file.
It can be used to obtained elastic moduli as a function of time to check their
convergence. If this option is used, -fgap/--frame-gap is an essential option.
Expand Down Expand Up @@ -100,7 +100,7 @@
"""Elastic Properties type.
Two keywords are accepted: "BST" or "ST".
* "BST" : Bending-Stretching-Twisting --- All motions are considered
" "ST" : Stretching-Twisting --- Bending motions are ignored.
* "ST" : Stretching-Twisting --- Bending motions are ignored.
WARNING: For accurate calculation of bending motions, DNA structures in trajectory must
be superimposed on a reference structure (See Publication's Method Section).
Expand All @@ -116,21 +116,19 @@
bpEndHelp=\
"""Last BP/BPS of DNA upto which parameter will be extracted.
If it is not given, parameter for only a single bp/s given with -bs/--bp-start
option will be extracted.
If it is not given, last basepair or base-step will be considered.
"""

frameGapHelp=\
"""Number of frames to skip for next calculation
"""Number of frames to skip for next calculation
When calculating elastic modulus as a function of time, this option will determine
the time-gap between each calculation point. Lower the time-gap, slower will be the
calculation.
"""

paxisHelp=\
"""Axis parallel to global helical-axis
"""Principal axis parallel to global helical-axis
Three keywords are accepted: "X", "Y" and "Z". Only require when bending motions are
included in the calculation.
Expand Down Expand Up @@ -171,10 +169,10 @@ def main():
showErrorAndExit(parser, "No Input File!!!\n")

# Determine file-extension type
fileType = 'text'
inputFileExtension = os.path.splitext(inputFile)[1]
if inputFileExtension in ['.h5', '.hdf5', 'hdf']:
fileType = 'hdf5'
if inputFileExtension not in ['.h5', '.hdf5', '.hdf']:
showErrorAndExit(parser, "File should be in HDF5 (h5, hdf5 or hdf file extension) format...\n")


# Total number of base-pair
firstBP = args.firstBP
Expand All @@ -184,42 +182,14 @@ def main():
else:
totalBP = args.totalBP

# Check for input parameter
if args.parameter is None:
showErrorAndExit(parser, "No Parameter name!!!\n")

if not ( (args.parameter in dnaMD.basePairTypeParameters) or \
(args.parameter in dnaMD.baseStepTypeParameters) ):
parser.print_help()
print("\n===== ERROR =======")
print('Unknown parameter name "{0}"!!!\n'.format(args.parameter))
print("Accepted parameters are as follows:")
count = 1
for parameter in dnaMD.basePairTypeParameters:
print('{0}. "{1}"'.format(count, parameter))
count += 1
for parameter in dnaMD.baseStepTypeParameters:
print('{0}. "{1}"'.format(count, parameter))
count += 1
print("\n===================")
sys.exit(-1)

# Output file
if args.outputFile is None:
showErrorAndExit(parser, "No output File!!!\n")

paramType = dnaMD.getParameterType(args.parameter)
toMinusBP = 1
if paramType == 'bps':
toMinusBP = 2

# Determine start and end-bp
toMinusBP = 2
startBP = args.startBP
if startBP is None:
startBP = firstBP
endBP = None
if args.endBP is not None:
endBP = args.endBP
endBP = args.endBP
if endBP is None:
endBP = firstBP + totalBP - toMinusBP

# Check consistency of start bp
if (startBP < firstBP) or (startBP > totalBP+firstBP-toMinusBP):
Expand All @@ -236,73 +206,117 @@ def main():
msg = 'The requested end bp {0} is out side of {1}-{2} range.'.format(endBP, firstBP, totalBP+firstBP-toMinusBP)
showErrorAndExit(parser, msg)

# Check if merge is required
if endBP is None:
merge = False
else:
merge = True
# Define DNA segement here
bp = [startBP, endBP]

# Check if merge-method is given
if merge and args.merge_method is None:
msg = 'No merging method is provided!!!!'
showErrorAndExit(parser, msg)
if args.esType == 'BST' and args.paxis is None:
showErrorAndExit(parser, 'To calculate bending, principal axis parallel to helical axis is required.')

# Store input file name
filename = None
if fileType == 'hdf5':
filename = inputFile
if args.vsTime is not None and args.frameGap is None:
showErrorAndExit(parser, 'Frame-gap is required when elasticity as a function of time is calculated.')

if endBP is not None:
bp = [startBP, endBP]
else:
bp = [startBP]
if args.err_type is not None and args.frameGap is None:
showErrorAndExit(parser, 'Frame-gap is required when error is calculated.')

# initialize DNA object
dna = dnaMD.DNA(totalBP, filename=filename, startBP=firstBP)
dna = dnaMD.dnaEY(totalBP, esType=args.esType, filename=inputFile, startBP=firstBP)

# Check if mask is in object
if dna.mask is not None:
if dna.dna.mask is not None:
masked = True
else:
masked = False

# Directly read the data from text file
if fileType == 'text':
# No support for smoothed axis, curvature and tangent
if args.parameter in ['helical x-axis smooth', 'helical y-axis smooth', 'helical z-axis smooth', 'helical axis curvature', 'helical axis tangent']:
print("\n===== ERROR =======")
print('Extraction of parameter "{0}" is only supported via input HDF5 file.'.format(args.parameter))
print("\n===== ERROR =======")
sys.exit(-1)
if args.vsTime is None:
out = ''
towrite = ''
mean_out = output_props(args.esType)
if args.esType == 'BST':
mean, result = dna.getStretchTwistBendModulus(bp, paxis=args.paxis, masked=True, matrix=False)
for i in range(4):
for j in range(4):
if j != 3:
out += '{0:>10.3f} '.format(result[i][j])
towrite += '{0:.3f},'.format(result[i][j])
else:
out += '{0:>10.3f}\n'.format(result[i][j])
towrite += '{0:.3f}\n'.format(result[i][j])

mean_out += '{0:>15.3f} '.format(mean[i])
mean_out += '\n'

if args.esType == 'ST':
mean, result = dna.getStretchTwistModulus(bp, masked=masked, matrix=False)
for i in range(2):
for j in range(2):
if j != 1:
out += '{0:.3f} '.format(result[i][j])
towrite += '{0:.3f},'.format(result[i][j])
else:
out += '{0:.3f}\n'.format(result[i][j])
towrite += '{0:.3f}\n'.format(result[i][j])

mean_out += '{0:>12.3f} '.format(mean[i])
mean_out += '\n'

print('=========== Elastic Modulus Matrix ===========\n')
print(out)
print('=========== ====================== ===========\n')

print('========================= Average Values ==========================\n')
print(mean_out)
print('=========== ====================== ====================== ===========')

if args.outputFile is not None:
with open(args.outputFile, 'w') as fout:
fout.write(towrite)

else:
time, modulus = dna.getModulusByTime(bp, args.frameGap, masked=masked, paxis=args.paxis, outFile=args.vsTime)

props_name = list(modulus.keys())
if args.err_type is not None:
error = dnaMD.get_error(time, list(modulus.values()), len(props_name), err_type=args.err_type, tool=args.tool)
sys.stdout.write("==============================================\n")
sys.stdout.write('{0:<16}{1:>14}{2:>14}\n'.format('Elasticity', 'Value', 'Error'))
sys.stdout.write("----------------------------------------------\n")
for i in range(len(props_name)):
sys.stdout.write('{0:<16}{1:>14.3f}{2:>14.3f}\n'.format(props_name[i], modulus[props_name[i]][-1], error[i]))
sys.stdout.write("==============================================\n\n")


def output_props(key):
out = ''
if key == 'BST':
props = ['Bending-1 Angle', 'Bending-2 Angle', 'Contour Length', 'Sum. Twist']
for i in range(4):
out += props[i] + ' '
out += '\n'

# Read and load parameter from file
dnaMD.setParametersFromFile(dna, inputFile, args.parameter, bp=bp)
if key == 'ST':
props = ['Contour Length', 'Sum. Twist']
for i in range(2):
out += props[i] + ' '
out += '\n'

# Extract the input parameter for input DNA/RNA segment
values, density = dna.parameter_distribution(args.parameter, bp, bins=30, merge=merge, merge_method=args.merge_method, masked=masked)
return out

# Write the extracted data in a text file
fout = open(args.outputFile, 'w')
fout.write('# "{0}" \t Density\n'.format(args.parameter))
for i in range(len(values)):
fout.write("{0:.6}\t{1:.6}\n".format(values[i], density[i]))
fout.close()

def parseArguments():
parser = argparse.ArgumentParser(prog='dnaMD histogram',
description=description,
formatter_class=argparse.RawTextHelpFormatter)
parser = argparse.ArgumentParser(prog='dnaMD globalElasticity',
description=description,
formatter_class=argparse.RawTextHelpFormatter)

parser.add_argument('-i', '--input', action='store',
type=argparse.FileType('rb'), metavar='L-BP_cdna.dat',
type=argparse.FileType('rb'), metavar='parameter.h5',
dest='inputFile', required=False, help=inputFileHelp)

parser.add_argument('-o', '--output', action='store',
type=str, metavar='output.csv',
dest='outputFile', help=outputFileHelp)

parser.add_argument('-ot', '--output-time', action='store',
type=str, metavar='elasicity_vs_time.dat',
type=str, metavar='elasicity_vs_time.csv',
dest='vsTime', help=vsTimeHelp)

parser.add_argument('-tbp', '--total-bp', action='store',
Expand All @@ -329,7 +343,7 @@ def parseArguments():
metavar='frames-to-skip',
help=frameGapHelp)

parser.add_argument('-paxis', '--principle-axis', action='store',
parser.add_argument('-paxis', '--principal-axis', action='store',
type=str, metavar='X', default=None,
choices=['X', 'Y', 'Z'],
dest='paxis', help=paxisHelp)
Expand All @@ -340,14 +354,14 @@ def parseArguments():
dest='err_type', help=errorHelp)

parser.add_argument('-gt', '--gromacs-tool', action='store',
type=str, metavar='g_analyze', default='g_analyze',
type=str, metavar='gmx analyze', default='gmx analyze',
dest='tool', help=toolsHelp)

parser.add_argument('-fbp', '--first-bp', action='store',
type=int, metavar='1', default=1,
dest='firstBP', help=bpFirstHelp)

idx = sys.argv.index("histogram") + 1
idx = sys.argv.index("globalElasticity") + 1
args = parser.parse_args(args=sys.argv[idx:])

return parser, args
Expand Down

0 comments on commit b43e14d

Please sign in to comment.