Stacking of two neighboring non-base-paired bases

[ksripath]

Dear Dr. Lu and X3DNA users,

I am relatively new to X3DNA, and am trying to use the x3dna_md.rb and extract_par.rb to analyze stacking and sliding of two neighboring bases in my RNA MD simulations. I've successfully used x3dna_md.rb on PDB files of each frame of one of my MD simulations, but I'm thinking this might not be the best way to analyze my data.  I believe that since the two bases I'm interesting in monitoring are not in base pairs themselves, they would be missed by find_pair. I was thinking that the analyze script might be a way to go, but if I'm not mistaken, analyze also takes the output of find_pair for its input.

Any help that can be provided would be much appreciated! Please also do let me know if there are threads that have already addressed my concern; I have looked for these but may have missed them.

[xiangjun]

Welcome to join the 3DNA-user community!

Please note that the x3dna_md.rb and extract_par.rb pair has been replaced by the x3dna_ensemble script, which consolidates ensemble-processing functionality in 3DNA v2.1 under one umbrella. Please update your 3DNA installation to the latest version.

Regarding your specific question, could you provide a reproducible example? A single frame from your MD simulations can be used to illustrate your point. Show step-by-step what you want to achieve, what's missing from 3DNA, and we can start from there.

Xiang-Jun

[ksripath]

Dear Dr. Lu,

Thank you for your response.

I want to analyze two bases in a ribozyme, G1 and U24, neither of which are in any base pairs. I want to track the sliding interaction of these two bases across the total length of my trajectory, so that I can correlate the sliding to various other distances and angles that I measure across the length of the trajectory as well. I can see how one can extract the slide parameter for a base pair step, but I don't understand if and how one can extract the slide parameter for two unpaired bases.

First I run find_pair as my first step, and you can see that neither G1 or U24 is listed in the output, which I have attached.

Next, if I use the 'extract' option for x3dna_ensemble, then if my understanding is correct, I would first need to run 'analyze' whose input I believe is from find_pair. I don't think this will work for me because my two bases of interest aren't listed inthe find_pair output.

I also see that there is an 'analyze' option for x3dna_ensemble. I am not sure if this will work for me, but I haven't found the correct syntax to enter my range of PDB files (I have generated one PDB file per frame of the trajectory that I wan't to analyze.). I've tried each of the following 3 commands:

> x3dna_ensemble analyze snaps/10000/10/*
> x3dna_ensemble -v analyze snaps/10000/10/*
> x3dna_ensemble -h analyze snaps/10000/10/*

where "snaps/10000/10/" is the directory that contains the first 10 pdb files from the first 10 frames of my trajectory.

And this is the message I get every time:

> Specify only one of: ensemble|models|pattern|list|one

Please let me know if you need further detail, Dr. Lu, and thank you so much again!

[xiangjun]

Thank for providing further details. However, for your example to be reproducible, you need to provide a PDB file -- that's the starting point of all the following steps.

Basically, you need to run find_pair first to get the pairing info in a file (as the one you attached); the file may need to be modified manually if some desired pair is missing, as is most likely in your case. Then one runs x3dna_ensemble analyze to get the various 3DNA parameters, followed by x3dna_ensemble extract to pick up the parameters one is interested in. The first step in analyzing an ensemble (MD or NMR), i.e, preparing a correct bp file, is identical to that for a single structure, and it is the most crucial.

In the coming 3DNA JoVE paper, there is a protocol to illustrate the whole procedure. In the meantime, you may want to try the examples distributed with 3DNA to get familiar on how to run x3dna_ensemble.

Xiang-Jun

[ksripath]

Dear Dr. Lu,

Thank you very much for the details on x3dna_ensemble; I'm just wondering how one would specify the pdb files of the ensemble? Or if one runs x3dna_ensemble analyze and x3dna_ensemble extract in the same directory as the pdb files to be analyzed, will x3dna_ensemble just run through all the pdb files available? I'm sorry for these questions; I've looked in the x3dna-v2.1/examples/ensemble/md directory, but haven't found a README file, and so am a bit confused about the syntax. Please don't hesitate to direct me to other threads or other sources if they exist and I haven't found them yet.

I'm sorry for not attaching the PDB file before; I only just now understood what you meant. I've attached it to this post.

I was wondering if you could provide more information on how one could manually modify the output from find_pair? As I said, the G1 and U24 bases are not in base pairs themselves, so I'm not sure if adding them to the find_pair output would help.

Thank you very much again, and sorry for more questions.

[xiangjun]

Thanks for providing a PDB file which helps clarify the issue.

The find_pair program is working as designed. For your RNA structure, it identifies two helical regions (see the attached figure below --- [Note added on 2013-01-28: image deleted upon request]):

Code: [Select]

pdb.00001
pdb.out
    2         # duplex
   30         # number of base-pairs
    1    1    # explicit bp numbering/hetero atoms
   31    8  0 #    1 | ....>-:..31_:[.RG]G-----C[.RC]:...8_:-<....  0.23  0.08 16.57  8.85 -3.78
   32    7  0 #    2 | ....>-:..32_:[.RC]C-----G[.RG]:...7_:-<....  0.50  0.23 11.02  8.87 -3.49
   33    6  0 #    3 | ....>-:..33_:[.RA]A-----U[.RU]:...6_:-<....  0.15  0.03  4.87  8.77 -4.55
   34    5  0 #    4 | ....>-:..34_:[.RA]A-----U[.RU]:...5_:-<....  0.14  0.09  6.89  8.85 -4.32
   35    4  0 #    5 | ....>-:..35_:[.RG]G-----C[.RC]:...4_:-<....  0.84  0.37  4.54  8.71 -3.19
   36    3  0 #    6 | ....>-:..36_:[.RC]C-----G[.RG]:...3_:-<....  0.84  0.79 15.85  8.75 -1.79
   37    2  0 #    7 | ....>-:..37_:[.RU]U-**--G[.RG]:...2_:-<....  3.01  0.85 17.27  8.99  2.57
   38   23  0 #    8 | ....>-:..38_:[.RG]G-----C[.RC]:..23_:-<....  0.45  0.30 12.77  8.96 -3.31
   39   22  0 #    9 | ....>-:..39_:[.RG]G-----C[.RC]:..22_:-<....  0.50  0.33 33.29  8.85 -2.19
   42   40  0 #   10 | ....>-:..42_:[.RA]A-**+-G[.RG]:..40_:-<....  7.23  1.99 47.35  7.57 10.57
   43   61  0 #   11 | ....>-:..43_:[.RA]A-**--G[.RG]:..61_:-<....  1.60  0.88 20.41 10.70  1.38
   44   60  0 #   12 | ....>-:..44_:[.RC]C-----G[.RG]:..60_:-<....  0.34  0.04 24.77  8.89 -3.35
   45   59  0 #   13 | ....>-:..45_:[.RA]A-----U[.RU]:..59_:-<....  0.13  0.06 17.54  8.79 -3.87
   46   58  0 #   14 | ....>-:..46_:[.RU]U-----A[.RA]:..58_:-<....  0.11  0.07 12.05  8.78 -4.15
   47   57  0 #   15 | ....>-:..47_:[.RU]U-----A[.RA]:..57_:-<....  0.36  0.34 10.79  8.72 -3.42
   48   55  0 #   16 | ....>-:..48_:[.RC]C-----G[.RG]:..55_:-<....  0.40  0.02 28.47  8.84 -3.15
   49   54  0 #   17 | ....>-:..49_:[.RC]C-----G[.RG]:..54_:-<....  0.53  0.38 22.85  8.85 -2.56
   50   53  9 #   18 x ....>-:..50_:[.RG]G-**--A[.RA]:..53_:-<....  9.64  0.99 16.65  8.51 11.45
   10   73  0 #   19 | ....>-:..10_:[RG5]g-----c[RC3]:..73_:-<....  0.37  0.20 12.89  8.98 -3.58
   11   72  0 #   20 | ....>-:..11_:[.RG]G-----C[.RC]:..72_:-<....  0.97  0.26 20.72  8.78 -2.48
   12   71  0 #   21 | ....>-:..12_:[.RU]U-----A[.RA]:..71_:-<....  0.17  0.04 22.04  9.00 -3.65
   13   70  0 #   22 | ....>-:..13_:[.RC]C-----G[.RG]:..70_:-<....  0.46  0.18 17.81  9.02 -3.30
   14   69  0 #   23 | ....>-:..14_:[.RC]C-----G[.RG]:..69_:-<....  0.52  0.46 10.86  9.01 -3.01
   15   68  0 #   24 | ....>-:..15_:[.RG]G-----C[.RC]:..68_:-<....  0.33  0.29 12.02  8.95 -3.50
   16   67  0 #   25 | ....>-:..16_:[.RC]C-----G[.RG]:..67_:-<....  0.36  0.30 11.63  8.96 -3.45
   17   66  0 #   26 | ....>-:..17_:[.RA]A-----U[.RU]:..66_:-<....  0.18  0.00  4.99  8.79 -4.57
   18   30  0 #   27 | ....>-:..18_:[.RG]G-----C[.RC]:..30_:-<....  0.35  0.07 24.12  8.82 -3.30
   19   29  0 #   28 | ....>-:..19_:[.RC]C-----G[.RG]:..29_:-<....  0.49  0.20 24.09  8.88 -2.91
   20   28  0 #   29 | ....>-:..20_:[.RC]C-----G[.RG]:..28_:-<....  0.73  0.69 13.00  9.01 -2.25
   21   26  0 #   30 | ....>-:..21_:[.RU]U-**+-G[.RG]:..26_:-<....  1.30  0.34 16.05  9.43 -0.22
##### Base-pair criteria used:     4.00     0.00    15.00     2.50    65.00     4.50     7.80 [ O N]
##### 5 non-Watson-Crick base-pairs, and 2 helices (0 isolated bps)
##### Helix #1 (18): 1 - 18  ***broken O3'[i] to P[i+1] linkage***
##### Helix #2 (12): 19 - 30  ***broken O3'[i] to P[i+1] linkage***
G1 and U24 are not paired according to the default criteria. Indeed, as shown below and as mentioned in your initial post, G1 and U24 are stacking instead of pairing.

I want to analyze two bases in a ribozyme, G1 and U24, neither of which are in any base pairs. I want to track the sliding interaction of these two bases across the total length of my trajectory, so that I can correlate the sliding to various other distances and angles that I measure across the length of the trajectory as well. I can see how one can extract the slide parameter for a base pair step, but I don't understand if and how one can extract the slide parameter for two unpaired bases.

If you insist on finding the relative geometry of G1 to U24, you could play the following trick:

Code: [Select]

find_pair -s pdb.00001 pdb.00001.datThe -s option means to output a list of all nucleotides in the given PDB file. Since you are interested in only G1 vs U24, you can manually edit the above output file pdb.00001.dat to have the following content:

Code: [Select]

pdb.00001
pdb.outs
    1      # single helix
    2      # number of bases
    1    1 # explicit bp numbering/hetero atoms
    1      #     1 ....>-:...1_:[RG5]g
   24      #    24 ....>-:..24_:[.RU]U
Running "analyze pdb.00001.dat", you get in file pdb.outs the following section:

Code: [Select]

****************************************************************************
Local base step parameters
    step       Shift     Slide      Rise      Tilt      Roll     Twist
   1  g/U      -1.95     -2.62      1.18    -29.34    146.33     78.47
****************************************************************************
Try a few frames (snapshots) along your MD simulation trajectories to decide for yourself if that make sense -- certainly this usage is beyond 3DNA's 'normal' application range.

Xiang-Jun

[ksripath]

Dear Dr. Lu,

Thank you so much for your response. I did think that perhaps my request would be outside x3DNA's normal usage, and that some parameters, like the "Rise" parameter are meaningless. It does seem to me, however, that the slide parameter might still have meaning, as you suggest.

Thank you very much again, and have a great day.

[ksripath]

Dear Dr. Lu,

Since I am conducting this analysis in preparation for a manuscript, would you mind removing the image of the whole structure that you posted, in case a competitor might recognize it? I may be being overly cautious, but I would prefer that to the alternative.

Thank you for your consideration, and have a great day.

[xiangjun]

Upon your request, I've removed "the image of the whole structure". I understand your concern, and good luck with your manuscript!

Xiang-Jun

[ksripath]

Thank you so much, Dr. Lu! I really appreciate it, as well as all your help!

Have a great day!

[ksripath]

Dear Dr. Liu,

I have a follow-up question on this discussion: I've been able to use the protocol that you suggested; however, I'm wondering if you have a suggestion as to how to optimize it for analyzing each frame of my trajectory. I have about 25,004 frames in my MD trajectory, and I would have to run find_pair -s for each of those frames, modify by hand each of those outputs and then run analyze on each of those outputs. Is this possible?

Thank you so much, and have a great day.

[xiangjun]

This is where programming (scripting) skills come into play. You do not have to do this manually for each of your 25k+ frames, once you have written a purpose-speicifc script to automate the process. You may find it helpful to have a look of the Ruby script x3dna_ensemble and corresponding files under folder lib/.

Alternatively, Curves+ or other tools may better suit your needs.

Xiang-Jun

[ksripath]

Dear Dr. Liu,

Thank you very much for your response. I'm trying to put find_pair in a for loop to write out separate outputs for each of my pdb files. find_pair seems to be working fine, but my outputs aren't being written. I was wondering if there was a bug that you knew of that wouldn't allow me to use find_pair in such a fashion. I realize that my questions may be outside the scope of this forum; I just wanted to bring it up in case this was an unknown bug with find_pair. Please do let me know if this is the case, and I very much appreciate all your help.

Thanks again, and have a great day.

[xiangjun]

Reviewing the thread, you would recall that find_pair is not used for your special case. The input file is prepared manually, and then fed to analyze to get the result. It is more likely something is wrong in your for-loop than a bug in 3DNA. If you understand the protocol as we went through, then try to write a for loop that repeats only once to confirm its correctness.

As always, a reproducible example would have made things clear.

Xiang-Jun