srnas / barnaba

Analyse Nucleic Acids Structure and Simulations with baRNAba
GNU General Public License v3.0
38 stars 15 forks source link

Calculating Pseudotorsion angle #45

Closed mandar5335 closed 2 years ago

mandar5335 commented 2 years ago

Hi, I am planning to calculate RNA pseudotorsion angles (eta, theta, eta', theta', eta'', and theta'' Ref: https://x3dna.org/highlights/pseudo-torsions-to-simplify-the-representation-of-dna-rna-backbone-conformation). If I define a section using corresponding atomnames in "nucleic.py" , is it sufficient? Or do I need to change some other parts of the code also? Thanks in advance.

Best, Mandar Kulkarni

sbottaro commented 2 years ago

Hello,

One option could be to define a new function in nucleic.py similar to get_bb_torsion_idx, i.e. something of this sort:

 def get_eta_theta_torsion_idx(self):
        # the first index is the residue number, the second corresponds to 
        # the six torsions and the third index is the atom involved
        idxs = np.zeros((len(idx_residues),6,4),dtype=int)
        rr = []
        for j,i in enumerate(idx_residues):
              # eta
              idxs[j,0] = [res.atom("P").index, res.atom("XX").index, res.atom("XX").index, res.atom("XX").index]
...

And another one in functions.py similar to backbone_angles_traj:

def eta_theta_angles_traj(traj):
    top = traj.topology
    # initialize nucleic class
    nn = nucleic.Nucleic(top)
    all_idx,rr =  nn.get_eta_theta_torsion(residues)
...

Let me know if this helps

Sandro

mandar5335 commented 2 years ago

Hi Sandro, Thanks for the suggestions and detailed reply. I have successfully added eta, theta, eta', and theta' prime.

I have modified nucleic.py, functions.py, definitions.py, and commandline.py. These files are attached here. barnaba_pseudotor.tar.gz

However, in case of eta, I am facing a small issue. First value of eta should be nan or absent, but, I am getting value as it could be considering last residue [-1]. at i=0. Could you please let me know what I am missing here?


      for j,i in enumerate(idx_residues):

            res = self.ok_residues[i]
            rr.append(self.rna_seq[i])
            # eta
            try:
                res_minus = self.ok_residues[i-1]
                res_plus = self.ok_residues[i+1]

                #print res,res_minus, res_minus.chain.index,res.chain.index
                idxs_tmp  = [res_minus.atom("C4'").index,res.atom("P").index,res.atom("C4'").index,res_plus.atom("P").index]

                if (i!=0 and (res_minus.chain.index == res.chain.index)):
                   idxs[j,0] = idxs_tmp
                else: pass

                if(res_plus.chain.index == res.chain.index):
                   #print "mm", idxs_tmp
                   idxs[j,0] = idxs_tmp
                else:
                    pass
            except:
                pass

Thanks a lot once again.

Best, Mandar Kulkarni

mandar5335 commented 2 years ago

Hi, The above issue is solved by merging two if conditions as:

                if (i!=0 and (res_minus.chain.index == res.chain.index) and (res_plus.chain.index == res.chain.index)):
                   idxs[j,0] = idxs_tmp
                else: pass

The modified files are attached here. pseudotor_files1.zip

If one wants to calculate RNA pseudotorsions , eta, theta, eta', and theta', these files should be copied to barnaba folder. Then, pseudotorsion can be calculated using the command line or in jupyter-notebook. command line:

barnaba TORSION --pseudo --pdb foo.pdb

Jupyter Notebook: definitions.pseudo_angles This will calculate four pseudotorsions: ['eta', 'theta', 'eta_prime', 'theta_prime']

Thanks again @sbottaro

Best, Mandar Kulkarni