2
$\begingroup$

I looking at a particular protein structure called 2LY4 accessible from RSCB PDB website. The corresponding fasta sequence for that structure is this:

>2LY4_1|Chain A|High mobility group protein B1|Homo sapiens (9606)
GKGDPKKPRGKMSSYAFFVQTCREEHKKKHPDASVNFSEFSKKCSERWKTMSAKEKGKFEDMAKADKARYEREMKTYIPPKGE
>2LY4_2|Chain B|Cellular tumor antigen p53|Homo sapiens (9606)
MEEPQSDPSVEPPLSQETFSDLWKLLPENNVLSPLPSQAMDDLMLSPDDIEQWFTEDPGPDEAPRMPEAAPPVAPAPAAPTPAAPAPAPSWPL

And the PDB format file can be downloaded here (1.7MB). The entire is to large to be pasted here.

What I want to do is to extract the subset of PDB format based on the subsequence in fasta above. Namely Chain A starting from 1st residue to 30th residue.

GKGDPKKPRGKMSSYAFFVQTCREEHKKKH

How can I do that in R or Python?

$\endgroup$
4
$\begingroup$

In general, if you simply want to extract that part of the PDB file, you could loop over it (it's plain text) and check the fields you're interested in:

with open('2ly4.pdb') as pdb:
    for line in pdb:
        if line[:4] == 'ATOM':
            chain = line[21]
            res_idx = int(line[23:26])
            if chain == 'A' and 1 <= res_idx <= 30:
                print(line[:-1])
        else:
            print(line[:-1]) # alternatively, skip with continue

However, this doesn't quite work with the file you've linked because it actually contains ten different conformations (models). This is better handled with a proper parser like prody:

import prody

model = prody.parsePDB('2ly4.pdb', model=1)
extract = model.select('chain A resnum 1to30')
prody.writePDB('extract.pdb', extract)

This will also allow you to do more complex atom selection.

$\endgroup$
3
  • $\begingroup$ Neat! The downside is that it removes all the headers including important ones, such as CRYST1. $\endgroup$
    – marcin
    Jun 7 at 13:27
  • $\begingroup$ If you remove model=1 all models are written (but with a small change: in MODEL 3 the charge -1 that in the original file is on atom 390 (OE1 in Glu 25) is moved to atom 391 OE2). $\endgroup$
    – marcin
    Jun 7 at 14:00
  • $\begingroup$ Right, parsing the file by hand gives you more control over the output, at the expense of reinventing the wheel. $\endgroup$ Jun 8 at 1:24
3
$\begingroup$

Here is an example how to do it using another Python library with another selection syntax.
Requires pip install gemmi.

import gemmi

st = gemmi.read_structure('2ly4.pdb')
subset = gemmi.parse_cid('/*/A/1-30').copy_structure_selection(st)
subset.write_minimal_pdb('output.pdb')
$\endgroup$
1
  • $\begingroup$ As @marcin suggested, it's probably better not to rely on ad hoc parser and use functionalities already provided by specialized libraries. $\endgroup$
    – saiden
    Jun 7 at 14:09
2
$\begingroup$

If you need to do it only once, for this one PDB file, it's easier to do it manually. When you open the file you can see it has 10 models from NMR. Do you want to extract a subset of one model or of all 10? In this file the residues are numbered according to the sequence. So you find the line where residue 31 begins in the first model:

ATOM    500  N   PRO A  31      18.844  -9.699  -4.620  1.00  0.00           N

and delete all lines until

ENDMDL

Then you do the same for other models, or you just remove other models, depending on what you need. You may also want to remove headers – from the beginning of the file until the CRYST1 record.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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