11
$\begingroup$

I'm trying to find a programmatic way to automatically extract the following information from a PDB file:

  1. RNA sequence
  2. Secondary structure restraints in bracket format, e.g. . (( . ( . ) . ))

Does software exist that can take a PDB file as input and generate these two pieces of information?

e.g. file. 3NDB_ba.pdb

The nucleotide sequence is:

gUCUCGUCCCGUGGGGCUCGGCGGUGGGGGAGCAUCUCCUGUAGGGGAGAUGUAACCCCCUUUACCUGCCGAACCCCGCCAGGCCCGGAAGGGAGCAACGGUAGGCAGGACGUCGGCGCUCACGGGGGUGCGGGAC.

And the secondary structure :

.(((.(..(((((((((.(((((..(((((.(((((((((....)))))))))..)))))....((((((.....(((.....(((....))).....)))..)))))).))))))).))))))).....).))).

Improve this question
$\endgroup$
3
  • $\begingroup$ I suggest you have a look at Biopython. They have a PDB module: biopython.org/wiki/The_Biopython_Structural_Bioinformatics_FAQ I'm not familiar with PDB, so I don't know if this does what you want. $\endgroup$
    – bli
    Jun 7, 2017 at 8:39
  • $\begingroup$ Hi Peter, thanks for posting your question, and welcome to Bioinformatics Stack Exchange. I have fixed up the grammar and spelling of your question to make it easier for other people to read and answer. I am not familiar with the PDB file format, so it would be great if you could add in a bit more context / story around your question. For example, where are these PDB files from (i.e. public, or private data source)? Is there a reason why taking transcript sequences from NCBI is not appropriate? $\endgroup$
    – gringer ♦
    Jun 7, 2017 at 11:01
  • $\begingroup$ Could you give an example of a specific PDB file we can test on? $\endgroup$
    – terdon ♦
    Jun 7, 2017 at 11:21

3 Answers 3

Reset to default
8
$\begingroup$

I suggest you take a look at rna-pdb-tools we do way more than you need! :-) The tools can get you a sequence, secondary structure and much more using various algorithms, and all is well documented http://rna-pdb-tools.readthedocs.io/en/latest/

To get sequence http://rna-pdb-tools.readthedocs.io/en/latest/main.html#get-sequence 

$ rna_pdb_tools.py --get_seq 5_solution_1.pdb
> 5_solution_1.pdb A:1-576
CAUCCGGUAUCCCAAGACAAUCUCGGGUUGGGUUGGGAAGUAUCAUGGCUAAUCACCAUGAUGCAAUCGGGUUGAACACUUAAUUGGGUUAAAACGGUGGGGGACGAUCCCGUAACAUCCGUCCUAACGGCGACAGACUGCACGGCCCUGCCUCAGGUGUGUCCAAUGAACAGUCGUUCCGAAAGGAAG

or you can get sequence and secondary structure via x3dna (http://x3dna.org/)

[mm] py3dna$ git:(master) ✗ ./rna_x3dna.py test_data/1xjr.pdb
test_data/1xjr.pdb
 >1xjr nts=47 [1xjr] -- secondary structure derived by DSSR
 gGAGUUCACCGAGGCCACGCGGAGUACGAUCGAGGGUACAGUGAAUU
 ..(((((((...((((.((((.....))..))..))).).)))))))

The problem is not trivial. I hope rna-pdb-tools has stuff that can solve this issue for you (http://rna-pdb-tools.readthedocs.io/en/latest/want.html#module-rna_pdb_tools.utils.rna_x3dna.rna_x3dna)

Improve this answer
$\endgroup$
8
  • 1
    $\begingroup$ @MarcinMagnus two things: first, it will fail on weird file names which is why it is generally bad practice to parse the output of ls. ls is designed to produce output for humans, not machines and can't deal with odd filenames gracefully (although recent versions of GNU ls have improved on that). Second, it's just inelegant since it isn't needed here. A simple for f in *txt will do the same thing and is far more robust since it can deal with any file name, even those containing newlines. $\endgroup$
    – terdon ♦
    Jun 7, 2017 at 19:59
  • 1
    $\begingroup$ Also, I apologize for my tone. The parsing ls thing is a classic meme in the *nix world and one of the more common pitfalls for novice shell programmers. I spend a lot of time in the *nix communities (I am a mod on the Unix & Linux and Ask Ubuntu sites) so I am probably more sensitive about this than I need to be. In any case, I didn't mean to disparage your work and my apologies if it came across that way. $\endgroup$
    – terdon ♦
    Jun 7, 2017 at 20:01
  • 2
    $\begingroup$ OK, I see it now ;-) I knew about abuse of cat, and echo. Now I see ls on the list porkmail.org/era/unix/award.html I didn't know that before. Thanks for pointing this out. I will fix the docs in all places ;-) thanks! I'm using Linux for more than 10y but still some things remind to be learned ;-) $\endgroup$
    – Marcin Magnus
    Jun 7, 2017 at 20:07
  • 1
    $\begingroup$ @MarcinMagnus thanks :) While you're at it, it would be great if you could also fix the missing quotes (for i in *; do command "$i"; done) since those are also essential when dealing with arbitrary file names. For more than you even wanted to know about the importance of quoting variables in *nix shells, see Security implications of forgetting to quote a variable in bash/POSIX shells and Why does my shell script choke on whitespace or other special characters?. $\endgroup$
    – terdon ♦
    Jun 7, 2017 at 20:12
  • 1
    $\begingroup$ I'll throw another suggestion into the mix regarding backticks, which can be replaced in sh scripts with $() (which is nestable), e.g. for i in $(ls *pdb) [but as per @terdon's comment the ls is redundant, so for i in *pdb would work better] $\endgroup$
    – gringer ♦
    Jun 7, 2017 at 23:39
4
$\begingroup$

Context

The PDB file format is a fixed-column file format designed in 1970s for storing structural models of macromolecules. The format has been around for long time, has many uses, and although it has official spec the files in circulation may not strictly conform to it. It always has a list of atoms with coordinates (the first two lines are added to stress that it's a fixed-column format, they are not part of the file):

         1         2         3         4         5         6         7         8
12345678901234567890123456789012345678901234567890123456789012345678901234567890
ATOM      1  N   VAL A   1       3.320  14.780   4.844  1.00 35.53           N
ATOM      2  CA  VAL A   1       3.577  16.239   4.984  1.00 35.39           C
ATOM      3  C   VAL A   1       4.896  16.398   5.727  1.00 30.43           C
ATOM      4  O   VAL A   1       5.143  15.702   6.732  1.00 30.51           O
ATOM      5  CB  VAL A   1       2.343  16.975   5.494  1.00 45.39           C
ATOM      6  CG1 VAL A   1       2.586  18.497   5.590  1.00 60.06           C
ATOM      7  CG2 VAL A   1       1.103  16.811   4.634  1.00 53.93           C
ATOM      8  N   LEU A   2       5.748  17.241   5.158  1.00 28.04           N
ATOM      9  CA  LEU A   2       7.116  17.471   5.661  1.00 24.31           C
ATOM     10  C   LEU A   2       7.166  18.490   6.792  1.00 24.00           C
...

or, to show some RNA:

ATOM     42  N3    G B   2       9.252  12.871  -1.168  1.00 36.98           N
ATOM     43  C4    G B   2       8.424  12.964  -2.233  1.00 34.09           C
ATOM     44  P     A B   3      10.376   8.321  -4.834  1.00 40.53           P
ATOM     45  OP1   A B   3      11.773   8.279  -5.364  1.00 38.97           O
ATOM     46  OP2   A B   3       9.396   7.283  -5.218  1.00 39.32           O
ATOM     47  O5'   A B   3      10.429   8.473  -3.211  1.00 36.55           O
ATOM     48  C5'   A B   3      11.698   8.232  -2.554  1.00 35.13           C

The Protein Data Bank -- international institution that archives and annotates structural models of biological molecules that anyone can deposit -- is now using mmCIF as the primary format. The fixed-column PDB format had inherent limitations (max. 99,999 atoms, but also it was hard to include additional info). They still generate PDB files if possible (i.e. except the largest structures), but the only format accepted now for depositions is mmCIF.

mmCIF has a rather not intuitive CIF syntax (think JSON or XML, but designed before XML), plus ontology defined by PDBx/mmCIF dictionary (think XML Schema). The content of mmCIF is divided into table-like categories with relations between tables (it was designed at the peak of RDMBS popularity) and is a bit harder to work with, so the old PDB format is more popular.

Sequence from PDB

The PDB format has a record called SEQRES that explicitly lists the sequence, for example (5NEO):

SEQRES   1 A   18    G   G   U   G   G   G   G   A   C   G   A   C   C          
SEQRES   2 A   18    C   C   A CBV   C

The mmCIF format has more details:

_entity_poly.pdbx_seq_one_letter_code       'GGUGGGGACGACCCCA(CBV)C' 
_entity_poly.pdbx_seq_one_letter_code_can   GGUGGGGACGACCCCACC 
_entity_poly.pdbx_strand_id                 A 

loop_
_entity_poly_seq.entity_id 
_entity_poly_seq.num 
_entity_poly_seq.mon_id 
_entity_poly_seq.hetero 
1 1  G   n 
1 2  G   n 
1 3  U   n 
1 4  G   n 
1 5  G   n 
1 6  G   n 
1 7  G   n 
1 8  A   n 
1 9  C   n 
1 10 G   n 
1 11 A   n 
1 12 C   n 
1 13 C   n 
1 14 C   n 
1 15 C   n 
1 16 A   n 
1 17 CBV n 
1 18 C   n

The sequence is conveniently extracted for you.

Additionally, _entity_poly_seq in mmCIF includes information about microheterogeneity, and the SEQRES record doesn't. It is relevant when the model has two alternative residues in the same place, because part of the sample had a point mutation.

PDB files that didn't come from PDB usually don't have the SEQRES record. You may get the sequence from the list of atoms (residue names in columns 18-20), but some residues may be missing in the atom list if atomic positions could not be determined.

So it should be easy to extract the sequence either way. For example, BioPython has module BioSeqIO with two PDB pseudo-formats:

  • pdb-seqres - Reads a Protein Data Bank (PDB) file to determine the complete protein sequence as it appears in the header (no dependencies).
  • pdb-atom - Uses Bio.PDB to determine the (partial) protein sequence as it appears in the structure based on the atom coordinate section of the file (requires NumPy for Bio.PDB).

I just tried BioPython 1.66 and while SeqIO.parse() can extract protein sequence, it fails with RNA, such as the 5NEO above. Ooops.

Secondary structure restraints

We may be thinking about different restraints here. I'll write about restraints used in refinement. The experimental data alone is normally not sufficient to refine a model, so one needs restraints that represent prior knowledge - lengths of atomic bonds, angles, planarity restraints, restraints based on local similarity to other structures that were determined from higher resolution data, etc, whatever can help to make a sensible model. Refinement programs (such as Refmac, BUSTER, Phenix.refine) try to fit the model to the data and satisfy geometrical restraints at the same time.

Macromolecular crystallographic software is quite fragmented and one normally uses several different programs in the process. It happens that the programs to generate secondary structure restraints are separate from the actual refinement programs. (Secondary structure restraints are less essential than restraints for covalent bonds). CCP4 has a program called LibG that makes DNA/RNA restraints for Refmac. Phenix has a program called phenix.secondary_structure_restraints, or you can use a server from UCSC.

These restraints don't have a "bracket format", but they explicitly specify expected distances and angles between atoms. For example (for phenix.refine):

  bond { 
    action = *add 
    atom_selection_1 = chain  A and resid    1  and name  O6 
    atom_selection_2 = chain  A and resid   18  and name  N4 
    distance_ideal = 2.91
    sigma = 0.1
  }

or (for Refmac):

exte dist first chain A resi 16 ins  . atom N6 second chain A resi 3 ins . atom O4 value 2.94 sigma 0.15 type 1 
exte dist first chain A resi 16 ins  . atom N1 second chain A resi 3 ins . atom N3 value 2.84 sigma 0.1 type 1 
exte torsion first chain A resi 16 ins  . atom C2 second chain A resi 16 ins . atom N1 third chain A resi 3 ins . atom N3 fourth chain A resi 3 ins . atom C4 value 180 sigma 15 type 1
Improve this answer
$\endgroup$
1
$\begingroup$

If you know how to use python:

1 - Download the ModeRNA module from here and install

2 - from a python IDLE execute:

from moderna import *
m = load_model('file.pdb', 'A') #A is the chain
seq = m.get_sequence()
sec = m.get_secstruc()

the variable sec store the secondary structure in dot-bracket notation.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.