<p>Here is an example how to do it using another Python library with another selection syntax.<br />
Requires <code>pip install gemmi</code>.</p>
<pre class="lang-py prettyprint-override"><code>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')
</code></pre>

CC BY-SA 4.0
2021-06-07T14:09:36.270Z
As @marcin suggested, it's probably better not to rely on ad hoc parser and use functionalities already provided by specialized libraries.
saiden
CC BY-SA 4.0
2021-06-07T14:04:53.660Z
2021-06-07T14:04:53.660Z
marcin
<p>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:</p>
<pre><code>ATOM    500  N   PRO A  31      18.844  -9.699  -4.620  1.00  0.00           N
</code></pre>
<p>and delete all lines until</p>
<pre><code>ENDMDL
</code></pre>
<p>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 <code>CRYST1</code> record.</p>

CC BY-SA 4.0
2021-06-07T06:42:26.870Z
2021-06-07T06:42:26.870Z
marcin
<p>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:</p>
<pre class="lang-py prettyprint-override"><code>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 &lt;= res_idx &lt;= 30:
                print(line[:-1])
        else:
            print(line[:-1]) # alternatively, skip with continue
</code></pre>
<p>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 <a href="http://prody.csb.pitt.edu/tutorials/structure_analysis/pdbfiles.html" rel="nofollow noreferrer">prody</a>:</p>
<pre class="lang-py prettyprint-override"><code>import prody

model = prody.parsePDB('2ly4.pdb', model=1)
extract = model.select('chain A resnum 1to30')
prody.writePDB('extract.pdb', extract)
</code></pre>
<p>This will also allow you to do more complex atom selection.</p>

CC BY-SA 4.0
2021-06-07T13:27:01.870Z
Neat! The downside is that it removes all the headers including important ones, such as `CRYST1`.
marcin
CC BY-SA 4.0
2021-06-07T14:00:27.587Z
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).
marcin
CC BY-SA 4.0
2021-06-08T01:24:08.990Z
Right, parsing the file by hand gives you more control over the output, at the expense of reinventing the wheel.
Jan Wilamowski
CC BY-SA 4.0
2021-06-07T07:07:04.823Z
2021-06-07T07:07:04.823Z
Jan Wilamowski
<p>I looking at a particular protein structure called <a href="https://www.rcsb.org/structure/2LY4" rel="nofollow noreferrer"><strong>2LY4</strong></a> accessible from <a href="https://www.rcsb.org/structure/2LY4" rel="nofollow noreferrer">RSCB PDB</a> website.
The corresponding fasta sequence for that structure is <a href="https://www.rcsb.org/fasta/entry/2LY4" rel="nofollow noreferrer">this</a>:</p>
<pre><code>&gt;2LY4_1|Chain A|High mobility group protein B1|Homo sapiens (9606)
GKGDPKKPRGKMSSYAFFVQTCREEHKKKHPDASVNFSEFSKKCSERWKTMSAKEKGKFEDMAKADKARYEREMKTYIPPKGE
&gt;2LY4_2|Chain B|Cellular tumor antigen p53|Homo sapiens (9606)
MEEPQSDPSVEPPLSQETFSDLWKLLPENNVLSPLPSQAMDDLMLSPDDIEQWFTEDPGPDEAPRMPEAAPPVAPAPAAPTPAAPAPAPSWPL
</code></pre>
<p>And the PDB format file can be downloaded <a href="https://files.rcsb.org/download/2LY4.pdb" rel="nofollow noreferrer">here</a> (1.7MB). The entire is to large to be pasted here.</p>
<p>What I want to do is to extract the subset of PDB format based on the subsequence in fasta above.
Namely <code>Chain A</code> starting  from 1st residue to 30th residue.</p>
<pre><code>GKGDPKKPRGKMSSYAFFVQTCREEHKKKH
</code></pre>
<p>How can I do that in R or Python?</p>

<r><python><fasta><protein-structure><pdb>
How to extract subset of protein structure (PDB format) file based on a subsequence of that protein
# How to extract subset of protein structure (PDB format) file based on a subsequence of that protein

(Asked by: username_0 on 2021-06-07)

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?

**Tags:** r, python, fasta, protein-structure, pdb

---

## 3 answers

### Answer by username_3 on 2021-06-07. Score: 4 | Accepted answer

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.

> **Comments:**
> **username_1**: Neat! The downside is that it removes all the headers including important ones, such as `CRYST1`.  
> **username_1**: 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).  
> **username_3**: Right, parsing the file by hand gives you more control over the output, at the expense of reinventing the wheel.  

---

### Answer by username_1 on 2021-06-07. Score: 3

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')

> **Comments:**
> **username_2**: As @marcin suggested, it's probably better not to rely on ad hoc parser and use functionalities already provided by specialized libraries.  

---

### Answer by username_1 on 2021-06-07. Score: 2

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.


---

