local helical axis vector for strongly bent DNA

[adihau]

Hello everyone,

I was analyzing a molecular dynamics simulation (http://www.amber.utah.edu/FF-DNA-bsc1-OL15/3GGI/opc/bsc0/ ) with do_x3dna and noticed, that in the HelAxis_g.dat file for most of the frames one line of data is missing. (sometimes even two lines are missing or replaced by '----'). The DNA contains 10 base-pairs, so there should be 9 helical axis vectors, since they are calculated as base-step parameters.

When I run single frames of the simulation with X3DNA, and look at the local helical axis vector there, I find four different cases.
1. "number of base-pairs: 10" ;  all helical axis vectors are calculated.
2. "number of base-pairs: 9"  ; the helical axis vector for one base-pair step is missing.
3. "number of base-pairs: 8"  ; the helical axis vectors for two base-pair steps are missing.
4. "This structure has broken O3'(i) to P(i+1) linkages"  and   "number of base-pairs: 9"  ; 
     the helical axis vector of one base-pair step is missing and in another line the data is replaced by '----'.

My explanation for this is, that if the program detects a broken base-pair, it refuses to calculate any vectors and prints '----' instead. Also it seems, that if the DNA is deformed or broken too much, then the program doesn't recognize this as a base-pair anymore and thinks the DNA is shorter than it actually is.

My question is now: Is there any way to force the program to still calculate a helical axis vector, even if the DNA is strongly bent and/ or broken?
And if yes, does maybe someone know, how to implement this into the do_x3dna program?

Thank you very much for your help,
Adina

[xiangjun]

Hi Adina,

Thanks for posting on the 3DNA Forum. Without seeing a concrete example, I can only guess based on experience. Let the structure (a single frame) be named sample.pdb, try the following to see if you can get what you want:

Code: [Select]

find_pair sample.pdb | analyze -c
Note the -c option which treats the detect double helix as a continous unit.

HTH,

Xiang-Jun

[adihau]

Hi Xiang-Jun,

thanks a lot for your fast reply! And sorry for not attaching the files.
I used the following command 

Code: [Select]

find_pair F29.pdb stdout | analyzeThe results, I get with that are F29-analyze.out and F30-analyze.out

The command you suggested, helps in the case of F29, where the '----' lines are calculated now (F29-analyze-c.out), but there is no change for F30.

The problem is, that this DNA has 10 base-pairs, but in F29-analyze.out it says, there are 9 base-pairs (And in F30-analyze.out it says there are 8 base-pairs). Why does X3DNA not recognize one (two) of the base-pairs?

Thanks for your help,
Adina

[xiangjun]

Hi Adina,

Thanks for your follow-up. Your cases occur frequently, especially in the field of MD simulations. Please have a look of the FAQ entry "How to fix missing (superfluous) base pairs identified by find_pair?".

Using your F30 as an example, the two terminal "base pairs" are too far off normal, as can be seen clearly from the attached image. Similar story happens for F29. For such cases, the 3DNA find_pair program won't detect these "presumed" base pairs. This is not a bug, but a feature by design. Users should pay close attention to these missing pairs.

You can force 3DNA to take them as pairs by manually editing the output from "find_pair" before feeding it into "analyze". The find_pair/analyze two-step procedure has been implemented in 3DNA specifically to allow for such flexibility.

HTH,

Xiang-Jun

[adihau]

Thank you very much, I was searching for something like this.

I have one more question.
Is there any possibility to find out, which of those criteria was the one, which caused, that the bases were not matched as pair? And can I find out what values those geometric quantities have?
Because then it would be easier to adjust the criteria in a way that all base-pairs are treated as such.

Thank you,
Adina

[xiangjun]

Check the source code to get to the bottom of the technical details.

For the analysis of NMR ensembles or MD trajectories, however, manually editing 'find_pair'-generated base-pair file may be the preferred way to go, as documented for the 'x3dna_ensemble' script distributed with 3DNA.

Xiang-Jun

[adihau]

Hi Xiang-Jun,

thank you for your answer. I still couldn't figure out a way, so that base-pairs are recognized, even if they are broken. Since I have a long trajectory, editing the 'find_pair'-generated base-pair file is maybe not the best option in this case.

I tried altering the criteria, but somehow there are still base-pairs not recognized, although I changed the criteria to quite high values. I will attach a few frames, with which I have problems. Do you know, what I did wrong? In those frames are broken base-pairs, but it should be possible to still match them, I think.

The criteria I changed were the following:
Section 3: for find_pair program
#   maximum distance between base origins  (was 15 originally)
<max_dorg>40.0</max_dorg>

#   maximum vertical base separation, i.e., |stagger|  (was 2.5 originally)
<max_dv>20</max_dv>

#   maximum angle between base normals (in range 0..90)  (was 65 originally)
<max_plane_angle>90.0</max_plane_angle>

#   minimum distance between RN9/YN1 base atoms   (was 4.5 originally)
<min_dNN>1</min_dNN>

#   distance criterion for helix break   (was 7.8 originally)
<helix_break>50</helix_break>


Thank you very much for your time and your help,
Adina

[xiangjun]

Imagine if you do not have access to 'find_pair', how would manually provide the base-pairing information from a PDB file to the 3DNA 'analyze' program? No changes in the criteria can fit every possible cases, especially for extremely deformed 'base pairs'. In your F29.pdb and F30.pdb files, for example, the two terminal 'pairs' may not be called base pairs at all. That's why I suggested you using the result from 'find_pair' as a starting point, and make changes manually as needed, and then feed the revised pairing info to 'analyze'.

Two of your new attached files (Frame25574_13.pdb and Frame46393_9.pdb) are more common, where one base-pair is not detected by 'find_pair' in its default settings. Using Frame25574_13.pdb as an example,

find_pair Frame25574_13.pdb

will give the following result:

Frame25574_13.pdb
Frame25574_13.out
    2         # duplex
   17         # number of base-pairs
    1     1    # explicit bp numbering/hetero atoms
    1    36   0 #    1 | ....>X:...1_:[DG5]G-----C[DC3]:..36_:X<....   1.11   0.58  20.02   8.77  -1.73
    2    35   0 #    2 | ....>X:...2_:[.DC]C-----G[.DG]:..35_:X<....   0.18   0.03  20.28   9.06  -3.74
    3    34   0 #    3 | ....>X:...3_:[.DA]A-----T[.DT]:..34_:X<....   0.34   0.29  13.33   9.12  -3.41
    4    33   0 #    4 | ....>X:...4_:[.DC]C-----G[.DG]:..33_:X<....   0.58   0.33   2.85   8.97  -3.62
    5    32   0 #    5 | ....>X:...5_:[.DG]G-----C[.DC]:..32_:X<....   0.43   0.15  11.47   9.00  -3.70
    7    30   0 #    6 | ....>X:...7_:[.DA]A-----T[.DT]:..30_:X<....   0.41   0.32  13.65   9.21  -3.26
    8    29   0 #    7 | ....>X:...8_:[.DC]C-----G[.DG]:..29_:X<....   0.55   0.48   8.48   8.98  -3.06
    9    28   0 #    8 | ....>X:...9_:[.DG]G-----C[.DC]:..28_:X<....   0.46   0.04  17.63   9.10  -3.58
   10    27   0 #    9 | ....>X:..10_:[.DA]A-----T[.DT]:..27_:X<....   0.39   0.04  10.30   9.03  -4.01
   11    26   0 #   10 | ....>X:..11_:[.DA]A-----T[.DT]:..26_:X<....   0.78   0.77  13.39   8.77  -2.02
   12    25   0 #   11 | ....>X:..12_:[.DC]C-----G[.DG]:..25_:X<....   0.81   0.80  12.58   9.12  -1.96
   13    24   0 #   12 | ....>X:..13_:[.DG]G-----C[.DC]:..24_:X<....   0.55   0.18  18.67   9.09  -3.15
   14    23   0 #   13 | ....>X:..14_:[.DA]A-----T[.DT]:..23_:X<....   0.52   0.28  17.64   8.98  -3.03
   15    22   0 #   14 | ....>X:..15_:[.DA]A-----T[.DT]:..22_:X<....   0.47   0.30  24.22   8.63  -2.72
   16    21   0 #   15 | ....>X:..16_:[.DC]C-----G[.DG]:..21_:X<....   0.52   0.27  10.14   9.00  -3.43
   17    20   0 #   16 | ....>X:..17_:[.DG]G-----C[.DC]:..20_:X<....   0.56   0.55   4.96   9.11  -3.09
   18    19   0 #   17 | ....>X:..18_:[DC3]C-----G[DG5]:..19_:X<....   0.33   0.19  16.06   9.22  -3.49
##### Base-pair criteria used:     4.00     0.00    15.00     2.50    65.00     4.50     7.80 [ O N]
##### 0 non-Watson-Crick base-pairs, and 1 helix (0 isolated bps)
##### Helix #1 (17): 1 - 17  ***broken O3' to P[i+1] linkage***

Note that the presumed base-pair between DA6 and DT31 is missed. See the the attached image for why it is not detected by 'find_pair'. If you want to have it the the output of 'analyze', simply add it in, as shown below. Note that the number of base-pairs has been increased from 17 to 18 to account for the additional pair.

Frame25574_13.pdb
Frame25574_13.out
    2         # duplex
   18         # number of base-pairs
    1     1    # explicit bp numbering/hetero atoms
    1    36   0 #    1 | ....>X:...1_:[DG5]G-----C[DC3]:..36_:X<....   1.11   0.58  20.02   8.77  -1.73
    2    35   0 #    2 | ....>X:...2_:[.DC]C-----G[.DG]:..35_:X<....   0.18   0.03  20.28   9.06  -3.74
    3    34   0 #    3 | ....>X:...3_:[.DA]A-----T[.DT]:..34_:X<....   0.34   0.29  13.33   9.12  -3.41
    4    33   0 #    4 | ....>X:...4_:[.DC]C-----G[.DG]:..33_:X<....   0.58   0.33   2.85   8.97  -3.62
    5    32   0 #    5 | ....>X:...5_:[.DG]G-----C[.DC]:..32_:X<....   0.43   0.15  11.47   9.00  -3.70
    6    31   0 
    7    30   0 #    6 | ....>X:...7_:[.DA]A-----T[.DT]:..30_:X<....   0.41   0.32  13.65   9.21  -3.26
    8    29   0 #    7 | ....>X:...8_:[.DC]C-----G[.DG]:..29_:X<....   0.55   0.48   8.48   8.98  -3.06
    9    28   0 #    8 | ....>X:...9_:[.DG]G-----C[.DC]:..28_:X<....   0.46   0.04  17.63   9.10  -3.58
   10    27   0 #    9 | ....>X:..10_:[.DA]A-----T[.DT]:..27_:X<....   0.39   0.04  10.30   9.03  -4.01
   11    26   0 #   10 | ....>X:..11_:[.DA]A-----T[.DT]:..26_:X<....   0.78   0.77  13.39   8.77  -2.02
   12    25   0 #   11 | ....>X:..12_:[.DC]C-----G[.DG]:..25_:X<....   0.81   0.80  12.58   9.12  -1.96
   13    24   0 #   12 | ....>X:..13_:[.DG]G-----C[.DC]:..24_:X<....   0.55   0.18  18.67   9.09  -3.15
   14    23   0 #   13 | ....>X:..14_:[.DA]A-----T[.DT]:..23_:X<....   0.52   0.28  17.64   8.98  -3.03
   15    22   0 #   14 | ....>X:..15_:[.DA]A-----T[.DT]:..22_:X<....   0.47   0.30  24.22   8.63  -2.72
   16    21   0 #   15 | ....>X:..16_:[.DC]C-----G[.DG]:..21_:X<....   0.52   0.27  10.14   9.00  -3.43
   17    20   0 #   16 | ....>X:..17_:[.DG]G-----C[.DC]:..20_:X<....   0.56   0.55   4.96   9.11  -3.09
   18    19   0 #   17 | ....>X:..18_:[DC3]C-----G[DG5]:..19_:X<....   0.33   0.19  16.06   9.22  -3.49
##### Base-pair criteria used:     4.00     0.00    15.00     2.50    65.00     4.50     7.80 [ O N]
##### 0 non-Watson-Crick base-pairs, and 1 helix (0 isolated bps)
##### Helix #1 (17): 1 - 17  ***broken O3' to P[i+1] linkage***

Does this solve your problem? If so, please repeat the above procedures for the other two cases (Frame46393_9.pdb and Frame72614_2.pdb), showing clearly which pair is missing using a graphic image as I did for the DA6--DT31 pair. This would improve your understanding of the concept, and benefit other users.

Xiang-Jun

[adihau]

Hello Xiang-Jun,

thank you very much for your detailed answer and your patience. This indeed helps me, and I understand that changing the criteria does not help in these cases.

So for everyone having the same problem, I used the following code:

Code: [Select]

find_pair myframe.pdb myframe.inpthis will create a file 'myframe.inp', which I modified in the way Xiang-Jun suggested in the answer above. Then

Code: [Select]

analyze myframe.inpcreates a file myframe.out with the results for all base-pairs including the broken one.
This worked in all cases so far, I will attach pictures of the problematic base-pairs of two more frames. I attached two pictures of the Frame72614_2 because it's not very clear to see, what is happening there in the first picture. The bases are not too far away from each other, but probably twisted too much.


Thank you again,
Adina

[xiangjun]

Hi Adina,

Thanks for your followup. Note that for your 'Frame72614_2', the terminal G1 and C36 are stacking instead of pairing.

Presumably, the three structures (Frame25574_13.pdb, Frame46393_9.pdb, and Frame72614_2.pdb) are along the trajectory of an MD simulation of the same DNA molecule. The base-pair information would be the same: you do not need to run 'find_pair' repeatedly for each frame. Simply run 'find_pair' once for a selected frame (e.g., the first one), and edit the resulting pairing info as needed. Then that pairing info can be used for the analysis of each and every frame along the trajectory.

Best regards,

Xiang-Jun

[Rajendra Kumar]

Hello,

I would like to point out that do_x3dna contains -ref option, which executes find_pair for only first frame, and same base-pairs information is used for whole trajectory. It is written in this link: http://do-x3dna.readthedocs.io/en/latest/do_x3dna_usage.html#ref.

With best regards,
Rajendra