Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-06-07T14:09:36.270Z
Answers.Text: As @marcin suggested, it's probably better not to rely on ad hoc parser and use functionalities already provided by specialized libraries.
Answers.UserDisplayName: saiden
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-06-07T14:04:53.660Z
Answers.LastActivityDate: 2021-06-07T14:04:53.660Z
Answers.UserDisplayName: marcin
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-06-07T06:42:26.870Z
Answers.LastActivityDate: 2021-06-07T06:42:26.870Z
Answers.UserDisplayName: marcin
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-06-07T13:27:01.870Z
Answers.Text: Neat! The downside is that it removes all the headers including important ones, such as `CRYST1`.
Answers.UserDisplayName: marcin
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-06-07T14:00:27.587Z
Answers.Text: 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).
Answers.UserDisplayName: marcin
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-06-08T01:24:08.990Z
Answers.Text: Right, parsing the file by hand gives you more control over the output, at the expense of reinventing the wheel.
Answers.UserDisplayName: Jan Wilamowski
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-06-07T07:07:04.823Z
Answers.LastActivityDate: 2021-06-07T07:07:04.823Z
Answers.UserDisplayName: Jan Wilamowski
Body: <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>
  
Tags: <r><python><fasta><protein-structure><pdb>
Title: How to extract subset of protein structure (PDB format) file based on a subsequence of that protein
ThreadText: # 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**](https://www.rcsb.org/structure/2LY4) accessible from [RSCB PDB](https://www.rcsb.org/structure/2LY4) website.
  The corresponding fasta sequence for that structure is [this](https://www.rcsb.org/fasta/entry/2LY4):
  
  ```
  >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](https://files.rcsb.org/download/2LY4.pdb) (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:
  
  ```py
  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](http://prody.csb.pitt.edu/tutorials/structure_analysis/pdbfiles.html):
  
  ```py
  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`.
  
  ```py
  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 @username_1 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.
  
  
  ---
  
  