Re: Ruby scripts for the analysis of MD simulation trajector

[temizna]

Dear Xiang-Jun,

Thank you for your scripts! I just tested the sample set and they seem to work fine. I am glad you coded an -all option. I personally like to run these scripts once and have all the information available at the same time. In addition, this is very similar to what Curves does as well. Next step would be to incorporate some distribution calculation.

I have one small comment. Instead of printing out model_1, model_2, etc.. We should just have the snapshot numbers, since an MD trajectory tends to have thousands of them.

I look forward to using these scripts in my analysis.

Cheers!
Alpay

[xiangjun]

Hi Alpay,

I just tested the sample set and they seem to work fine.

Thanks for being the first user to verify that the scripts are working as expected. I'd been expecting feedback from shahabshariati after I revised the scripts to v0.2 following his/her 'bug' report.

I added the [mono:1dzrze8o]-all[/mono:1dzrze8o] option simply because it makes sense to me. As mentioned previously, I have no MD experience. Nevertheless, I am glad to know that this functionality coincides with Curves.

Regarding the next step, what do you mean exactly "to incorporate some distribution calculation"? Some simple statistics, maybe?

Regarding your comment/suggestion,

Instead of printing out model_1, model_2, etc.. We should just have the snapshot numbers, since an MD trajectory tends to have thousands of them.

Do you mean to change the 1st column of each output parameter file from e.g., '[mono:1dzrze8o][red:1dzrze8o]model_10[/red:1dzrze8o][/mono:1dzrze8o]' to '[mono:1dzrze8o][red:1dzrze8o]10[/red:1dzrze8o][/mono:1dzrze8o]'? If that's the case, you can simply change line [mono:1dzrze8o]257[/mono:1dzrze8o] in function [mono:1dzrze8o]write_out_parameters()[/mono:1dzrze8o] of [mono:1dzrze8o]x3dna_md.rb[/mono:1dzrze8o] from
[pre:1dzrze8o]c0 = "model_#{entries[idx]}"[/pre:1dzrze8o]to [pre:1dzrze8o]c0 = entries[idx][/pre:1dzrze8o]
Thanks for your feedback. Please clarify so I can improve the scripts.

Xiang-Jun

[temizna]

Dear Xiang-Jun,

After I submitted my reply, I already changed the line to have the snapshots as numbers only in your script. Still, I think it is better to have the numbers as default. May be in your next version?

By distribution calculation I meant in addition to the means, standard deviations you already supply, it would be nice to have BI/BII distributions averaged over the trajectory, histograms of each parameter (per base pair/base pair step and over the whole structure).  I know there is no perfect analysis script/program out there, but the more the initial analysis does the better for the end user. It also encourages people to use it more!

One more suggestion: In the simulations there are always end effects. A general practice among MD people is to remove the first and last base pairs from the analysis to reduce these effects in the subsequent analysis. It would be great to have this option available as a choice in the script.

Cheers!
Alpay

[xiangjun]

Hi Alpay,

As always, thanks for your feedback and suggestions. Now specifically,

After I submitted my reply, I already changed the line to have the snapshots as numbers only in your script.  Still, I think it is better to have the numbers as default. May be in your next version?

Suggestion taken -- now in v0.3, the first column is just the snapshot number, without the [mono:2688raft]model_[/mono:2688raft] prefix.

By distribution calculation I meant in addition to the means, standard deviations you already supply, it would be nice to have BI/BII distributions averaged over the trajectory, histograms of each parameter (per base pair/base pair step and over the whole structure). I know there is no perfect analysis script/program out there, but the more the initial analysis does the better for the end user. It also encourages people to use it more!

Just to clarify, the scripts currently do not provide mean/std etc statistics -- they are intended to be calculated by the users, using R/Matlab/Octave/Excel  etc. Regarding BI/BII classification of backbone, histograms of parameters etc, I am putting them in my to-do list.

One more suggestion: In the simulations there are always end effects. A general practice among MD people is to remove the first and last base pairs from the analysis to reduce these effects in the subsequent analysis. It would be great to have this option available as a choice in the script.

Suggestion taken -- now v0.3 has a new [mono:2688raft][red:2688raft]-e[/red:2688raft][/mono:2688raft] option to accomplish just what you asked. Moreover, the two ends can be asymmetrical, e.g., [mono:2688raft][red:2688raft]-e 1 2[/red:2688raft][/mono:2688raft] to remove the first and the last two bps from the extracted parameter list.

Download v0.3 to have a try, and report back any issues you experience; and of course, any new suggestions!

Xiang-Jun

[temizna]

Thanks for the changes Xiang-Jun. I just ran the test, it looks good. Now I will move on to do a real (and long) analysis run of a cruciform. I will let you know how it turns out.
Regarding BI/BII, for my most recent work ( I am preparing the manuscript now) I used R to do most of the analysis after I generated the parameters with curves+. I calculated BI/BII myself using a variety of e-z cut-offs. I know most people use 0, but there was a 93 paper by Hartmann and Lavery which suggested 70. When I looked at the distribution of e-z from my simulations I saw that 50 was a better separation. I ended up using 0 to be consistent with previous publications. I just wanted to let you know how BI/BII turned out for me.

Alpay

[temizna]

Thanks for the reply Xiang-Jun.

I am trying to process a big trajectory file. 50000 models. When I run the script, I get the following error:
rb x3dna_md.rb -b ../cf/run1/cf.dat -e ../cf/run1/cf.1.tr.pdb
        ../cf/run1/cf.1.tr.pdb: with model numbers <= 0
x3dna_md.rb:96:in `block (2 levels) in parse_base_pair_parameters': undefined method `[]' for nil:NilClass (NoMethodError)
        from x3dna_md.rb:96:in `collect'
        from x3dna_md.rb:96:in `block in parse_base_pair_parameters'
        from x3dna_md.rb:95:in `each'
        from x3dna_md.rb:95:in `each_with_index'
        from x3dna_md.rb:95:in `parse_base_pair_parameters'
        from x3dna_md.rb:203:in `block in parse_3dna_output'
        from x3dna_md.rb:200:in `open'
        from x3dna_md.rb:200:in `parse_3dna_output'
        from x3dna_md.rb:281:in `block in process_ensemble_models'
        from x3dna_md.rb:272:in `each'
        from x3dna_md.rb:272:in `process_ensemble_models'
        from x3dna_md.rb:79:in `main'
        from x3dna_md.rb:497:in `<main>'
Process model #0 / 50000
-------------------------

It seems the script is getting the correct number of models as evidenced from the presence of a model_list.dat

I will try shorter trajectories to see if the same error exists.
Alpay

[xiangjun]

From the error message, I guess you have some setting problems, e.g., mismatch between base-pair file and snapshots, or the snapshot contains unknown nucleotide that needs to be specified in file '[mono:2fyzp058]baselist.dat[/mono:2fyzp058]' etc. Thus, it is likely that no sensible 3DNA output parameters file was generated.

It is always a good idea to try a simple example with only one or two snapshots to verify that everything is set up correctly; afterwards, it is only a timing issue between 100 or one million snapshots.

Depending on user feedbacks, I am planning to refine the scripts to check common pitfalls and output sensible error messages accordingly. The next v0.4 should be much more robust for real world applications.

Xiang-Jun

[temizna]

The problem is temp_model.out file. I removed everything from the tar file but the scripts, the bpfile.dat and the sample_md.pdb, then ran the script. It gave me an error. Then I ran analyze -c temp_model.inp temp_model.pdb and created temp_model.out
Running the ensemble option again, processed the whole trajectory. I tested this on my own trajectory with 50 snapshots and now thing seem to work.


HTH,
Alpay

[xiangjun]

Hi Alpay,

I have updated the Ruby scripts to v0.4 which fixed all the glitches that came to my mind.  Please have a try and report back how it goes.

Xiang-Jun

[temizna]

Hi Xiang-Jun

After my initial problems, the scripts are working nicely. I have only one problem/request. I noticed that the x3dna_md script is dumping the results after finishing processing all the snapshots/models. This means that it is keeping a lot in the memory. MD simulations tend to have thousands of snaphsots that needs to be processed. I can of course skip snapshots and process for example every 10, but I really do not want to throw away any information. Is it possible for the script to append to the output file after processing each model?

Thanks
Alpay

[xiangjun]

Hi Alpay,

Glad to hear your verification that "the scripts are working nicely". I am well aware of the memory issue you raised. By design, the scripts currently process each and every snapshot of an MD simulation and collect all the parameters along the way. As to how many snapshots can be handled, the answer would obviously depend on computer memory and the size (i.e. no. of bps) of a snapshot. I would imagine the scripts should work well for an MD simulation with 50k snapshots,12 bps each, in a computer with say 4GB memory. As an estimation, the output parameters file for the sample 12 bps/20 snapshots takes 76,926 bytes. With 50k snapshots, the space taken would be  (76,926 / 20) * 50000 = 192,300,000 bytes, which is ~183 MB.

HTH,

Xiang-Jun

[jorgeRNA]

hello friends

  I have a problem when I try run the Ruby scripts for the analysis of MD simulation trajectory

 I follow the instructions of README text.

EXAMPLES

        Change to the sub-directory 'examples/' under your installation of the
        scripts (e.g., x3dna_md_v0.6/examples), then:

                x3dna_md.rb -b bpfile.dat -e sample_md0.pdb -o my.out
          or
                ../x3dna_md.rb -b bpfile.dat -e sample_md0.pdb -o my.out
          and
                diff x3dna_md.out my.out

---------------------------------------------------------------------------------------------------------

~/3dna/x3dna-v2.1/bin/x3dna_md_v0.7/examples$ x3dna_md.rb -b bpfile.dat -e sample_md0.pdb -o my.out
x3dna_md.rb: not order found

--------------------------------------------------------------------------------------------------------
but, I don't know exactly what happen ...

also try the follow, and the messege is :

:~/3dna/x3dna-v2.1/bin$ ./x3dna_md.rb
X3DNA environmental variable not set

please, any help will be appreciated

jorge

 

[xiangjun]

Hi Jorge,

You have not downloaded or installed 3DNA v2.1 correctly. The script is now called x3dna_ensemble. You seems not that familiar with the installation procedure. Please ask a local computer expert for help to get started.

Or, post back in details what you have done including the version of 3DNA you download, its location, your OS. Did you run x3dna_setup? What's the output? The point is reproducibility -- we must see exactly where your problems are to help you out.

Xiang-Jun

[jorgeRNA]

ok,

so, the thing is that... is just that this script is now called x3dna_ensemble.

I had installed the x3dna correctly, just look this lines in my root (terminal)

~/3dna/x3dna-v2.1/bin$ ls
  anyhelix      comb_str   find_platform  o1p_o2p   regular_dna  step_hel              x3dna_setup
  blocview      ex_str     frame_mol      pdb2img   rotate_mol   x3dna_ensemble        x3dna_utils
alc2img   bp_order.dat  fiber      get_part       r3d_atom  stack2img    x3dna_md_v0.7
analyze   cehs          find_pair  mutate_bases   rebuild   std_base     x3dna_md_v0.7.tar.gz

entering into the directory

~/3dna/x3dna-v2.1/bin/x3dna_md_v0.7$ ls
examples  extract_par.rb  README  trollop.rb  x3dna_md.rb

but you tell me that it isn't the correct script and now it has other name

so now I will  check your help and report to the forum

thanks

jorge

[xiangjun]

Well, for verification, I have just installed 3DNA v2.1 both Mac OS X and Linux 64-bit on my computer. The content of the bin/ directory is:

Code: [Select]

bin [522] ls
alc2img*        cehs*           find_pair*      mutate_bases*   rebuild*        std_base*       x3dna_utils*
analyze*        comb_str*       find_platform*  o1p_o2p*        regular_dna*    step_hel*
anyhelix*       ex_str*         frame_mol*      pdb2img*        rotate_mol*     x3dna_ensemble*
blocview*       fiber*          get_part*       r3d_atom*       stack2img*      x3dna_setup*

What OS do you have? Please install 3DNA v2.1 in a new directory, and make sure you have the same content as above. If so, then go to $X3DNA/examples/ensemble/md/ and follow README file there.

Xiang-Jun