How to properly use x3dna_ensemble?

[oskolkov]

Dear Xiang-Jun and 3DNA users,

I am doing MD simulations of DNA segment using GROMACS 4.0.5 which produces trajectories as xtc and trr files, please find enclosed the initial structure pdb file, I did not manage to attach xtc trajectory (21 Mb) due to size restrictions. I want to analyze the trajectories in order to obtain twist, roll, tilt etc profiles. I found out that x3dna_ensemble script included into x3dna_2.1beta version can provide this information. Unfortunately, I have difficulty running this script and would like to ask you what exactly I should type in my bash in order to analyze my xtc/trr trajectories? For example, I tried

$  x3dna_ensemble convert trajectory.xtc
{:help=>false, :package=>"amber"}

or

$  x3dna_ensemble convert trajectory.xtc -p gromacs
{:help=>false, :package_given=>true, :package=>"gromacs"}

in order to convert MD trajectories to MODEL/ENDMDL format, which can be later analyzed by "x3dna_ensemble analyze", but without success. Do I need to somehow split my xtc/trr trajectory into separate pdb files and apply analyze afterwards? I am completely lost, please help me. Could you please write exact commands I should type in my bash? Thank you.

With best regards,
Nikolay.

**********************************************************************
Dr. Nikolay Oskolkov
Department of Clinical Sciences, Diabetes and Endocrinology,
CRC, Entr 72, building 91, level 12, University Hospital Malmo
SE-205 02 Malmo, Sweden

Email: Nikolay.Oskolkov@med.lu.se
Tel.:  +46 761463349

[xiangjun]

Hi Nikolay,

Thanks for using 3DNA, especially trying the new Ruby x3dna_ensemble utility. Sorry to hear that you are "completely lost". Regarding the convert sub-command, you are ahead of time:

Code: [Select]

x3dna_ensemble convert -h
...to be added for Amber/Gromacs/CHARMM etc
  --package, -p <s>:   Name of MD simulation package (default: amber)
         --help, -h:   Show this message
i.e., this functionality is currently not implemented yet. However, I'd be interested in adding such a converter for GROMACS if you provide me an example (shortened) trajectory file, and let me know the detailed description of the xtc/trr file formats. On the other hand, I believe (or guess) GROMACS should have a facility to convert its native trajectory file to the standard PDB MODEL/ENDMDL format.

Please let me know.

Xiang-Jun

[oskolkov]

Dear Xiang-Jun,

thank you very much for your reply! Yes, you are right, I checked and found out that Gromacs can actually devide an xtc-trajectory file into separate pdb files in MODEL/ENDML format. So in fact I can apply "analyze" command to each pdb file/trajectory and find mechanical properties (local ones) of genome segment at each moment of time. In order to get tilt, twist, roll etc. for the whole period of simulation (whole xtc trajectory), I have to read and average the output information from a number of out-files, since as far as I understand X3DNA can provide mechanical information only for one pdb-file but not for a group of pdb-files which form a trjectory, do I understand it correctly? Curves+ can do it for the whole trajectory, but is much more tricky to use compared to X3DNA

Regards,
Nikolay.

[xiangjun]

Hi Nikolay,

I am glad to hear that Gromacs provides the facility to output an ensemble of multiple models delineated by MODEL/ENDMDL. As of 3DNA v2.1beta (currently distributed), the x3dna_ensemble Ruby script can handle multiple structures in a PDB MODEL/ENDMDL ensemble:

Code: [Select]

x3dna_ensemble -h
------------------------------------------------------------------------                                                       
Utilities for the analysis and visualization of an ensemble
    Usage: x3dna_ensemble [-h|-v] sub-command [-h] [options]
    where sub-command must be one of:
        analyze -- analyze MODEL/ENDMDL delineated ensemble (NMR or MD)
        block_image -- generate a base block schematic image
        extract -- extract structural parameters after running 'analyze'
        reorient -- reorient models to a particular frame/orientation
------------------------------------------------------------------------
  --version, -v:   Print version and exit
     --help, -h:   Show this message

Check the $X3DNA/examples/ensemble/ directory for examples, and report back if you have any problem. Note that I still need to add a detailed documentation of the new x3dna_ensemble utility. However, it should be straightforward to play with, and I am always quicker in responding to user's questions than writing doc ...

Xiang-Jun

[oskolkov]

Hi Xiang-Jun,

everything seems to work perfectly. I used Gromacs' function trjconv to divide the xtc trjectory into separate snapshots, put them to a subdirectory, generated bpfile.dat using find_pair function of 3DNA and the initial pdb-conformation of my DNA segment, and applied x3dna_ensemble analyze script. The output file bp_step.par contains all the local inter- (twist, roll, tilt etc) and intra-base pair step parameters (shear, stretch etc.). As far as I understand, this is the result of averaging with respect to all snapshots, i.e. time averaged, please correct me if I am wrong.

Cheers!
Nikolay.

[xiangjun]

Hi Nikolay,

Glad to hear that you've made progress. However, from my understanding of what you described, it seems something is still not quite right. The fixed-name file bp_step.par contains only the parameters for a single structure (snapshot), not the whole ensemble. It is a bit more difficult to explain the details in text, so I would suggest you repeat the examples in the directory $X3DNA/examples/ensemble/md/:

Code: [Select]

cd $X3DNA/examples/ensemble/md/
x3dna_ensemble analyze -h
x3dna_ensemble extract -hOnce you understand how the examples work, you should be able to apply the same idea to the analysis of your MD trajectories. Of course, if you have any questions, please do not hesitate to post back at the forum.

Xiang-Jun

PS: command-line helpThe help page for x3dna_ensemble analyze

Code: [Select]

------------------------------------------------------------------------
Analyze a MODEL/ENDMDL delineated ensemble of NMR structures or MD
trajectories. All models must correspond to different conformations of
the same molecule. A template base-pair input file, generated with
'find_pair' and corrected manually as necessary, must be provided.

Usage:
        x3dna_ensemble analyze options
Examples:
        x3dna_ensemble analyze -b bpfile.dat -e sample_md0.pdb
             # 21 models (0-20); output (default): 'ensemble_example.out'
             # also generate 'model_list.dat', see example below
        x3dna_ensemble analyze -b bpfile.dat -m model_list.dat -o ensemble_example2.out
             # diff ensemble_example.out ensemble_example2.out

        x3dna_ensemble analyze -b bpfile.dat -p 'pdbdir/model_*.pdb' -o ensemble_example3.out
             # note to quote the -p option; 20 models (1-20)
             # also generate 'pdb_list.dat', see example below
        x3dna_ensemble analyze -b bpfile.dat -l pdb_list.dat -o ensemble_example4.out
             # diff ensemble_example3.out ensemble_example4.out
             # note the order of the models: 1, 10..19, 2, 20, 3..9
Options:
------------------------------------------------------------------------
    --bpfile, -b <s>:   Name of file containing base-pairing info
   --outfile, -o <s>:   Output file (default: ensemble_example.out)
  --ensemble, -e <s>:   Ensemble delineated with MODEL/ENDMDL pairs
    --models, -m <s>:   File containing an explicit list of model numbers
   --pattern, -p <s>:   Pattern of model files to process (e.g., *.pdb)
      --list, -l <s>:   File containing an explicit list of models
          --info, -i:   Show only model info in the ensemble [with -e]
          --help, -h:   Show this message

The help page for x3dna_ensemble extract

Code: [Select]

------------------------------------------------------------------------
Extract 3DNA structural parameters of an ensemble of NMR structures or
MD trajectories, after running 'x3dna_ensemble analyze'. The extracted
parameters are intended to be exported into Excel, Matlab and R etc for
further data analysis/visualization.

Usage:
        x3dna_ensemble extract options
Examples:
        x3dna_ensemble extract -l
             # to see a list of all parameters
        x3dna_ensemble extract -p prop
             # for propeller, no need to specify full: -p pr suffices
             # -p 36 also fine (see above); use 'ensemble_example.out'
        x3dna_ensemble extract -p slide -s , -f ensemble_example3.out
             # comma separated, from file 'ensemble_example3.out'
        x3dna_ensemble extract -p roll -s ' ' -n -o roll.dat
             # space separated, no row-label, to file 'roll.dat'
        x3dna_ensemble extract -e 1 -p chi1
             # extract the chi torsion angle of strand I, but exclude
             # those from the two terminal base pairs. For comparison,
             # run also: x3dna_ensemble extract -p chi1
        x3dna_ensemble extract -a
             # extract all parameters, each in a separate file
Options:
------------------------------------------------------------------------
  --separator, -s <s>:   Separator for fields [\t] (default: )
   --par-name, -p <s>:   Name of parameter to extract
   --fromfile, -f <s>:   Parameters file (default: ensemble_example.out)
    --outfile, -o <s>:   File of selected parameter (default: stdout)
   --end-bps, -e <i+>:   Number of end pairs to ignore (default: 0, 0)
            --all, -a:   Extract all parameters into separate files
          --clean, -c:   Clean up parameter files by the -a option
           --list, -l:   List all parameters
        --no-1col, -n:   Delete the first (label) column
           --help, -h:   Show this message


[oskolkov]

Hi Xiang-Jun,

thank you for correcting me. Let me decribe step-by-step the commands I am using now

1)cd $X3DNA/examples/ensemble/md/

2) find_pair sample_md0.pdb bpfile.dat

which gives me bpfile the template base pair input file bpfile.dat. My bpfile.dat looks slightly different from the bpfile.dat provided in the examples/ensemble/md. My bpfile.dat has 13 base-pairs, not 12 as in examples/ensemble/md.

3) x3dna_ensemble analyze -b bpfile.dat -p 'pdbdir/model_*.pdb' -o ensemble.out
Process PDB file pdbdir/model_1.pdb     1 / 20
      5 [main] ruby 5696 child_info_fork::abort: address space needed by 'etc.so' (0x250000) is already occupied
      7 [main] ruby 5740 child_info_fork::abort: address space needed by 'etc.so' (0x250000) is already occupied
      7 [main] ruby 3340 child_info_fork::abort: address space needed by 'etc.so' (0x250000) is already occupied
Error in running command: 'C:\Packages\X3DNA/bin/find_pair -s temp_model.pdb nts_list.txt 2> msgfile'
Rerun the command w/o redirection for details

Very strange, now it gives me an error when I apply x3dna_ensemble analyze. It worked perfectly first time.

Regards,
Nikolay.

3)

[xiangjun]

Hi Nikolay,

The first two steps are fine. The third step is weird -- are you using MinGW or Cygwin version? This is a problem not necessarily specific to 3DNA, but obviously I'd like to see a solution to it.

Xiang-Jun

[oskolkov]

Yes, I am currently using Cygwin version

[oskolkov]

Dear Xiang-Jun,

I reinstalled X3DNA and this time it worked out. I got the following

$  x3dna_ensemble analyze -b bpfile.dat -p 'pdbdir/model_*.pdb' -o ensemble_example3.out
Process PDB file pdbdir/model_1.pdb     1 / 20
      5 [main] ruby 5276 child_info_fork::abort: address space needed by 'etc.so' (0x370000) is already occupied
Process PDB file pdbdir/model_10.pdb    2 / 20
Process PDB file pdbdir/model_11.pdb    3 / 20
      5 [main] ruby 4496 child_info_fork::abort: address space needed by 'etc.so' (0x370000) is already occupied
Process PDB file pdbdir/model_12.pdb    4 / 20
Process PDB file pdbdir/model_13.pdb    5 / 20
      5 [main] ruby 2316 child_info_fork::abort: address space needed by 'etc.so' (0x370000) is already occupied
      5 [main] ruby 5512 child_info_fork::abort: address space needed by 'etc.so' (0x370000) is already occupied
Process PDB file pdbdir/model_14.pdb    6 / 20
      7 [main] ruby 6132 child_info_fork::abort: address space needed by 'etc.so' (0x370000) is already occupied
Process PDB file pdbdir/model_15.pdb    7 / 20
      7 [main] ruby 6120 child_info_fork::abort: address space needed by 'etc.so' (0x370000) is already occupied
Process PDB file pdbdir/model_16.pdb    8 / 20
Process PDB file pdbdir/model_17.pdb    9 / 20
      5 [main] ruby 4668 child_info_fork::abort: address space needed by 'etc.so' (0x370000) is already occupied
Process PDB file pdbdir/model_18.pdb    10 / 20
      5 [main] ruby 4428 child_info_fork::abort: address space needed by 'etc.so' (0x370000) is already occupied
Process PDB file pdbdir/model_19.pdb    11 / 20
Process PDB file pdbdir/model_2.pdb     12 / 20
Process PDB file pdbdir/model_20.pdb    13 / 20
Process PDB file pdbdir/model_3.pdb     14 / 20
Process PDB file pdbdir/model_4.pdb     15 / 20
Process PDB file pdbdir/model_5.pdb     16 / 20
Process PDB file pdbdir/model_6.pdb     17 / 20
      5 [main] ruby 4240 child_info_fork::abort: address space needed by 'etc.so' (0x370000) is already occupied
      5 [main] ruby 4976 child_info_fork::abort: address space needed by 'etc.so' (0x370000) is already occupied
      5 [main] ruby 5384 child_info_fork::abort: address space needed by 'etc.so' (0x370000) is already occupied
    235 [main] ruby 4256 child_info_fork::abort: address space needed by 'etc.so' (0x370000) is already occupied
      6 [main] ruby 4876 child_info_fork::abort: address space needed by 'etc.so' (0x370000) is already occupied
      5 [main] ruby 4128 child_info_fork::abort: address space needed by 'etc.so' (0x370000) is already occupied
      7 [main] ruby 5096 child_info_fork::abort: address space needed by 'etc.so' (0x370000) is already occupied
      7 [main] ruby 2872 child_info_fork::abort: address space needed by 'etc.so' (0x370000) is already occupied
      5 [main] ruby 2896 child_info_fork::abort: address space needed by 'etc.so' (0x370000) is already occupied
Process PDB file pdbdir/model_7.pdb     18 / 20
      5 [main] ruby 2392 child_info_fork::abort: address space needed by 'etc.so' (0x370000) is already occupied
Process PDB file pdbdir/model_8.pdb     19 / 20
      6 [main] ruby 6040 child_info_fork::abort: address space needed by 'etc.so' (0x370000) is already occupied
Process PDB file pdbdir/model_9.pdb     20 / 20



Then I typed
$ x3dna_ensemble extract -a , -f ensemble_example3.out

and got several .out files with base-pair step parameters.

[xiangjun]

Glad to know that you've made some progress. However, the message as shown below still bothers me:

ruby 5276 child_info_fork::abort: address space needed by 'etc.so' (0x370000) is already occupied

A quick Google search turns out quite a few hits concerning Ruby, and Cygwin on Windows. Will following Cygwin FAQ #4.44 "How do I fix fork() failures?" solve your problem? Specifically, the following sentence seems to address this problem:

Read the 'rebase' package README in /usr/share/doc/rebase/, and follow the instructions there to run 'rebaseall'

Additionally, what version of Ruby are you using?

Code: [Select]

ruby -v
I am switching from Perl to Ruby as the scripting language for 3DNA, hopefully this won't cause practical issues for Windows users.

Xiang-Jun

[oskolkov]

Hi Xiang-Jun,

thank you for the link, I will try to fix the problem with Cygwin. The version of Ruby I am using is

ruby 1.8.7 (2012-02-08 patchlevel 358) [i386-cygwin]

Nikolay

[oskolkov]

Hi Xiang-Jun,

finally, I figured out that the main parameter I am intereested in doing my DNA simulation is the bending persistence length, which is related to the DNA bending angle theta as <theta^2>=2L/L_p, where L is the contour length of the DNA segment. Thus, I have to calculate theta^2 averaged over all my MD trajectories.

Here come my question. Doing "analyze" of a single snapshot, I get an output file where (at the end) there is an information about helical axis vectors

Position (Px, Py, Pz) and local helical axis vector (Hx, Hy, Hz)
         for each dinucleotide step

BUT, when I do "x3dna_ensamble analyze" for multiple snapshots, I do not see helical axis vectors in the output file. Is that true, am I doing everything correctly? Is "x3dna_ensamble analyze" still not able to provide the helical axis vectors information despite "analyze" can do it? Looking forward to your reply.

Cheers!
Nikolay.
 

[xiangjun]

BUT, when I do "x3dna_ensamble analyze" for multiple snapshots, I do not see helical axis vectors in the output file.

The information you need, "Position (Px, Py, Pz) and local helical axis vector (Hx, Hy, Hz) for each dinucleotide step",  is not parsed in the current version of x3dna_ensemble. I will get it added soon and keep you updated.

Xiang-Jun
 

[oskolkov]

Thank you very much Xiang-Jun! That would be really helpful. Looking forward to the updated version.

[xiangjun]

Hi Nikolay,

Check the updated version dated 2012apr16 on the download page, and let me know if you have any problem.

Xiang-Jun

[oskolkov]

Hi Xiang-Jun,

thank you very much for this update! I have tested it, it worked perfectly. Such an enormous relief! Thank you!

Best wishes,
Nikolay.